home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / ARGONET / PD / MATHS / RLAB / RLAB125.ZIP / !RLaB / test < prev    next >
Text File  |  1995-11-18  |  73KB  |  3,194 lines
  1. //
  2. // Test Suite for RLaB.
  3. // This test suite is supposed to test as much as possible, and still
  4. // be runable on most platforms. This means we cannot do graphics or
  5. // other platform specific stuff like piping. However, that is OK,
  6. // since we are mostly interested in assuring the builder that RLaB
  7. // will produce correct numerical answers.
  8.  
  9. // We use randsvd, which in turn uses RLaB's random number
  10. // generator. We try to be carefull and seed the random number
  11. // generator so that each user will get similar results (or at least
  12. // similar inputs to the tests).
  13.  
  14. // Test the other style comments 
  15. 1 + 2; // A simple statement with trailing comment.
  16. #  Optional RLaB comment style.
  17. 1 + 2; #  A simple statement with trailing comment.
  18. %  Optional Matlab comment style.
  19. 1 + 2; %  A simple statement with trailing comment.
  20.  
  21. srand (SEED = 10);   // Seems to produce reasonable results
  22. rand("default");
  23.  
  24. tic();    // Start timing the tests...
  25.  
  26. //
  27. // Test Parameters and some functions we will need later
  28. //
  29.  
  30. pi = 4.0*atan(1.0);
  31. X = 3;    // should be 3 (heuristic).
  32.  
  33. //
  34. // Compute machine epsilon
  35. //
  36.  
  37. epsilon = function() 
  38. {
  39.   eps = 1.0;
  40.   while((eps + 1.0) != 1.0) 
  41.   {
  42.     eps = eps/2.0;
  43.   }
  44.   return 2*eps;
  45. };
  46.  
  47. eps = epsilon();
  48.  
  49. eye = function( m , n ) 
  50. {
  51.   if (!exist (n))
  52.   {
  53.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  54.     new = zeros (m[1], m[2]);
  55.     N = min ([m[1], m[2]]);
  56.   else
  57.     if (class (m) == "string" || class (n) == "string") {
  58.       error ("eye(), string arguments not allowed");
  59.     }
  60.     if (max (size (m)) == 1 && max (size (n)) == 1)
  61.     {
  62.       new = zeros (m[1], n[1]);
  63.       N = min ([m[1], n[1]]);
  64.     else
  65.       error ("matrix arguments to eye() must be 1x1");
  66.     }
  67.   }
  68.   for(i in 1:N)
  69.   {
  70.     new[i;i] = 1.0;
  71.   }
  72.   return new;
  73. };
  74.  
  75. symm = function( A )
  76. {
  77.   return (A + A')./2;
  78. };
  79.  
  80. //-------------------------------------------------------------------//
  81.  
  82. // Synopsis:    Pascal matrix.
  83.  
  84. // Syntax:    P = pascal ( N )
  85.  
  86. // Description:
  87.  
  88. //    The Pascal matrix of order N: a symmetric positive definite
  89. //    matrix with integer entries taken from Pascal's triangle. The
  90. //    Pascal matrix is totally positive and its inverse has integer
  91. //    entries.  Its eigenvalues occur in reciprocal pairs. COND(P)
  92. //    is approximately 16^N/(N*PI) for large N. PASCAL(N,1) is the
  93. //    lower triangular Cholesky factor (up to signs of columns) of
  94. //    the Pascal matrix.   It is involutary (is its own
  95. //    inverse). PASCAL(N,2) is a transposed and permuted version of
  96. //    PASCAL(N,1) which is a cube root of the identity.
  97.  
  98. //      References:
  99. //      R. Brawer and M. Pirovino, The linear algebra of the Pascal matrix,
  100. //           Linear Algebra and Appl., 174 (1992), pp. 13-23 (this paper
  101. //           gives a factorization of L = PASCAL(N,1) and a formula for the
  102. //           elements of L^k).
  103. //      S. Karlin, Total Positivity, Volume 1, Stanford University Press,
  104. //           1968.  (Page 137: shows i+j-1 choose j is TP (i,j=0,1,...).
  105. //                   PASCAL(N) is a submatrix of this matrix.)
  106. //      M. Newman and J. Todd, The evaluation of matrix inversion programs,
  107. //           J. Soc. Indust. Appl. Math., 6(4):466--476, 1958.
  108. //      H. Rutishauser, On test matrices, Programmation en Mathematiques
  109. //           Numeriques, Editions Centre Nat. Recherche Sci., Paris, 165,
  110. //           1966, pp. 349-365.  (Gives an integral formula for the
  111. //           elements of PASCAL(N).)
  112. //      J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  113. //           Birkhauser, Basel, and Academic Press, New York, 1977, p. 172.
  114. //      H.W. Turnbull, The Theory of Determinants, Matrices, and Invariants,
  115. //           Blackie, London and Glasgow, 1929.  (PASCAL(N,2) on page 332.)
  116.  
  117. //    This file is a translation of pascal.m from version 2.0 of
  118. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  119. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  120.  
  121. //-------------------------------------------------------------------//
  122.  
  123. pascal = function ( n , k )
  124. {
  125.   local(n, k)
  126.  
  127.   if (!exist (k)) { k = 0; }
  128.  
  129.   P = diag( (-1).^[0:n-1] );
  130.   P[; 1] = ones(n,1);
  131.  
  132.   //  Generate the Pascal Cholesky factor (up to signs).
  133.  
  134.   for (j in 2:n-1)
  135.   {
  136.     for (i in j+1:n)
  137.     {
  138.       P[i;j] = P[i-1;j] - P[i-1;j-1];
  139.     }
  140.   }
  141.  
  142.   if (k == 0)
  143.   {
  144.     P = P*P';
  145.   else if (k == 2) {
  146.     P = rot90(P,3);
  147.     if (n/2 == round(n/2)) { P = -P; }
  148.   }}
  149.  
  150.   return P;
  151. };
  152.  
  153. //-------------------------------------------------------------------//
  154.  
  155. // Synopsis:    Random, orthogonal upper Hessenberg matrix.
  156.  
  157. // Syntax:      H = ohess ( N )
  158.  
  159. // Description:
  160.  
  161. //      H is an N-by-N real, random, orthogonal upper Hessenberg
  162. //      matrix. Alternatively, H = OHESS(X), where X is an arbitrary
  163. //      real N-vector (N > 1) constructs H non-randomly using the
  164. //      elements of X as parameters. In both cases H is constructed
  165. //      via a product of N-1 Givens rotations. 
  166.  
  167. //      Note: See Gragg (1986) for how to represent an N-by-N
  168. //      (complex) unitary Hessenberg matrix with positive subdiagonal
  169. //      elements in terms of 2N-1 real parameters (the Schur
  170. //      parametrization). This M-file handles the real case only and
  171. //      is intended simply as a convenient way to generate random or
  172. //      non-random orthogonal Hessenberg matrices.
  173.  
  174. //      Reference:
  175. //      W.B. Gragg, The QR algorithm for unitary Hessenberg matrices,
  176. //      J. Comp. Appl. Math., 16 (1986), pp. 1-8.
  177.  
  178. //    This file is a translation of ohess.m from version 2.0 of
  179. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  180. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  181.  
  182. //-------------------------------------------------------------------//
  183.  
  184. ohess = function ( x )
  185. {
  186.   local (x)
  187.   global (pi)
  188.  
  189.   if (any (imag (x))) { error("Parameter must be real."); }
  190.  
  191.   n = max(size(x));
  192.  
  193.   if (n == 1)
  194.   {
  195.     //  Handle scalar x.
  196.     n = x;
  197.     x = rand(n-1, 1)*2*pi;
  198.     H = eye(n,n);
  199.     H[n;n] = sign(rand());
  200.   else
  201.     H = eye(n,n);
  202.     H[n;n] = sign(x[n]) + (x[n]==0);   // Second term ensures H[n;n] nonzero.
  203.   }
  204.  
  205.   for (i in n:2:-1)
  206.   {
  207.     // Apply Givens rotation through angle x(i-1).
  208.     theta = x[i-1];
  209.     c = cos(theta);
  210.     s = sin(theta);
  211.     H[ [i-1, i] ;] = [ c*H[i-1;]+s*H[i;] ;
  212.                        -s*H[i-1;]+c*H[i;] ];
  213.   }
  214.  
  215.   return H;
  216. };
  217.  
  218. //-------------------------------------------------------------------//
  219.  
  220. // Synopsis:    Random matrix with pre-assigned singular values.
  221.  
  222. // Syntax:    R = randsvd (N, KAPPA, MODE, KL, KU)
  223.  
  224. // Description:
  225.  
  226. //    R is a (banded) random matrix of order N with COND(A) = KAPPA
  227. //    and singular values from the distribution MODE.
  228.  
  229. //      N may be a 2-vector, in which case the matrix is N(1)-by-N(2).
  230. //      Available types:
  231. //             MODE = 1:   one large singular value,
  232. //             MODE = 2:   one small singular value,
  233. //             MODE = 3:   geometrically distributed singular values,
  234. //             MODE = 4:   arithmetically distributed singular values,
  235. //             MODE = 5:   random singular values with unif. dist. logarithm.
  236.  
  237. //      If omitted, MODE defaults to 3, and KAPPA defaults to SQRT(1/EPS).
  238. //      If MODE < 0 then the effect is as for ABS(MODE) except that in the
  239. //      original matrix of singular values the order of the diagonal entries
  240. //      is reversed: small to large instead of large to small.
  241. //      KL and KU are the lower and upper bandwidths respectively; if they
  242. //      are omitted a full matrix is produced.
  243. //      If only KL is present, KU defaults to KL.
  244. //      Special case: if KAPPA < 0 then a random full symmetric positive
  245. //                    definite matrix is produced with COND(A) = -KAPPA and
  246. //                    eigenvalues distributed according to MODE.
  247. //                    KL and KU, if present, are ignored.
  248.  
  249. //      This routine is similar to the more comprehensive Fortran routine xLATMS
  250. //      in the following reference:
  251. //      J.W. Demmel and A. McKenney, A test matrix generation suite,
  252. //      LAPACK Working Note #9, Courant Institute of Mathematical Sciences,
  253. //      New York, 1989.
  254.  
  255. //    This file is a translation of randsvd.m from version 2.0 of
  256. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  257. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  258.  
  259. // Dependencies
  260. #   rfile bandred
  261. #   rfile qmult
  262.  
  263. //-------------------------------------------------------------------//
  264.  
  265. randsvd = function (n, kappa, mode, kl, ku)
  266. {
  267.   local (n, kappa, mode, kl, ku)
  268.  
  269.   if (!exist (kappa)) { kappa = sqrt(1/eps); }
  270.   if (!exist (mode))  { mode = 3; }
  271.   if (!exist (kl))    { kl = n-1; }    // Full matrix.
  272.   if (!exist (ku))    { ku = kl; }     // Same upper and lower bandwidths.
  273.  
  274.   if (abs(kappa) < 1) {
  275.     error("Condition number must be at least 1!");
  276.   }
  277.  
  278.   posdef = 0;
  279.   if (kappa < 0)            // Special case.
  280.   {
  281.     posdef = 1;
  282.     kappa = -kappa;
  283.   }
  284.  
  285.   p = min(n);
  286.   m = n[1];        // Parameter n specifies dimension: m-by-n.
  287.   n = n[max(size(n))];
  288.  
  289.   if (p == 1)        // Handle case where A is a vector.
  290.   {
  291.     rand("normal", -10, 10);
  292.     A = rand(m, n);
  293.     A = A/norm(A, "2");
  294.     return A;
  295.   }
  296.  
  297.   j = abs(mode);
  298.  
  299.   // Set up vector sigma of singular values.
  300.   if (j == 3)
  301.   {
  302.     factor = kappa^(-1/(p-1));
  303.     sigma = factor.^[0:p-1];
  304.  
  305.   else if (j == 4) {
  306.     sigma = ones(p,1) - (0:p-1)'/(p-1)*(1-1/kappa);
  307.  
  308.   else if (j == 5) {    // In this case cond(A) <= kappa.
  309.     rand("uniform", 0, 1)
  310.     sigma = exp( -rand(p,1)*log(kappa) );
  311.  
  312.   else if (j == 2) {
  313.     sigma = ones(p,1);
  314.     sigma[p] = 1/kappa;
  315.  
  316.   else if (j == 1) {
  317.     sigma = ones(p,1)./kappa;
  318.     sigma[1] = 1;
  319.  
  320.   }}}}}
  321.  
  322.  
  323.   // Convert to diagonal matrix of singular values.
  324.   if (mode < 0) {
  325.     sigma = sigma[p:1:-1];
  326.   }
  327.  
  328.   sigma = diag(sigma);
  329.  
  330.   if (posdef)        // Handle special case.
  331.   {
  332.     Q = qmult(p);
  333.     A = Q'*sigma*Q;
  334.     A = (A + A')/2;    // Ensure matrix is symmetric.
  335.     return A;
  336.   }
  337.  
  338.   if (m != n) 
  339.   {
  340.     sigma[m; n] = 0;    // Expand to m-by-n diagonal matrix.
  341.   }
  342.  
  343.   if (kl == 0 && ku == 0) // Diagonal matrix requested - nothing more to do.
  344.   {
  345.     A = sigma;
  346.     return A;
  347.   }
  348.  
  349.   // A = U*sigma*V, where U, V are random orthogonal matrices from the
  350.   // Haar distribution.
  351.   A = qmult(sigma');
  352.   A = qmult(A');
  353.  
  354.   if (kl < n-1 || ku < n-1)    // Bandwidth reduction.
  355.   {
  356.    A = bandred(A, kl, ku);
  357.   }
  358.  
  359.   rand("default");
  360.   return A;
  361. };
  362.  
  363. //-------------------------------------------------------------------//
  364.  
  365. // Synopsis:    Band reduction by two-sided unitary transformations.
  366.  
  367. // Syntax:    bandred ( A , KL, KU )
  368.  
  369. // Description:
  370.  
  371. //    bandred(A, KL, KU) is a matrix unitarily equivalent to A with
  372. //    lower bandwidth KL and upper bandwidth KU (i.e. B(i,j) = 0 if
  373. //    i > j+KL or j > i+KU).  The reduction is performed using
  374. //    Householder transformations. If KU is omitted it defaults to
  375. //    KL. 
  376.  
  377. //    Called by randsvd.
  378. //    This is a `standard' reduction.  Cf. reduction to bidiagonal
  379. //    form prior to computing the SVD.  This code is a little
  380. //    wasteful in that it computes certain elements which are
  381. //    immediately set to zero! 
  382. //
  383. //      Reference:
  384. //      G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  385. //      Johns Hopkins University Press, Baltimore, Maryland, 1989.
  386. //      Section 5.4.3.
  387.  
  388. //    This file is a translation of bandred.m from version 2.0 of
  389. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  390. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  391.  
  392. // Dependencies
  393. #   rfile house
  394.  
  395. //-------------------------------------------------------------------//
  396.  
  397. bandred = function ( A , kl , ku )
  398. {
  399.   local (A, kl, ku)
  400.  
  401.   if (!exist (ku)) { ku = kl; else ku = ku; }
  402.  
  403.   if (kl == 0 && ku == 0) {
  404.     error ("You''ve asked for a diagonal matrix.  In that case use the SVD!");
  405.   }
  406.  
  407.   // Check for special case where order of left/right transformations matters.
  408.   // Easiest approach is to work on the transpose, flipping back at the end.
  409.  
  410.   flip = 0;
  411.   if (ku == 0)
  412.   {
  413.     A = A';
  414.     temp = kl; kl = ku; ku = temp; flip = 1;
  415.   }
  416.  
  417.   m = A.nr; n = A.nc; 
  418.  
  419.   for (j in 1 : min( min(m, n), max(m-kl-1, n-ku-1) ))
  420.   {
  421.     if (j+kl+1 <= m)
  422.     {
  423.        ltmp = house(A[j+kl:m;j]);
  424.        beta = ltmp.beta; v = ltmp.v;
  425.        temp = A[j+kl:m;j:n];
  426.        A[j+kl:m;j:n] = temp - beta*v*(v'*temp);
  427.        A[j+kl+1:m;j] = zeros(m-j-kl,1);
  428.     }
  429.  
  430.     if (j+ku+1 <= n)
  431.     {
  432.        ltmp = house(A[j;j+ku:n]');
  433.        beta = ltmp.beta; v = ltmp.v;
  434.        temp = A[j:m;j+ku:n];
  435.        A[j:m;j+ku:n] = temp - beta*(temp*v)*v';
  436.        A[j;j+ku+1:n] = zeros(1,n-j-ku);
  437.     }
  438.   }
  439.  
  440.   if (flip) {
  441.     A = A';
  442.   }
  443.  
  444.   return A;
  445. };
  446.  
  447. //-------------------------------------------------------------------//
  448.  
  449. // Synopsis:    Lehmer matrix - symmetric positive definite.
  450.  
  451. // Syntax:      A = lehmer ( N )
  452.  
  453. // Description:
  454.  
  455. //      A is the symmetric positive definite N-by-N matrix with
  456. //                     A(i,j) = i/j for j >= i.
  457. //      A is totally nonnegative.  INV(A) is tridiagonal, and explicit
  458. //      formulas are known for its entries. 
  459.  
  460. //      N <= COND(A) <= 4*N*N.
  461.  
  462. //      References:
  463. //        M. Newman and J. Todd, The evaluation of matrix inversion
  464. //           programs, J. Soc. Indust. Appl. Math., 6 (1958), pp. 466-476.
  465. //        Solutions to problem E710 (proposed by D.H. Lehmer): The inverse
  466. //           of a matrix, Amer. Math. Monthly, 53 (1946), pp. 534-535.
  467. //        J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  468. //           Birkhauser, Basel, and Academic Press, New York, 1977, p. 154.
  469.  
  470. //    This file is a translation of lehmer.m from version 2.0 of
  471. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  472. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  473.  
  474. //-------------------------------------------------------------------//
  475.  
  476. lehmer = function ( n )
  477. {
  478.   local (n)
  479.   global (tril)
  480.  
  481.   A = ones(n,1)*(1:n);
  482.   A = A./A';
  483.   A = tril(A) + tril(A,-1)';
  484.  
  485.   return A;
  486. };
  487.  
  488. //-------------------------------------------------------------------//
  489.  
  490. // Synopsis:    Pre-multiply by random orthogonal matrix.
  491.  
  492. // Syntax:    B = qmult ( A )
  493.  
  494. // Description:
  495.  
  496. //    B is Q*A where Q is a random real orthogonal matrix from the
  497. //    Haar distribution, of dimension the number of rows in
  498. //    A. Special case: if A is a scalar then QMULT(A) is the same as
  499.  
  500. //        qmult(eye(a))
  501.  
  502. //       Called by RANDSVD.
  503.  
  504. //       Reference:
  505. //       G.W. Stewart, The efficient generation of random
  506. //       orthogonal matrices with an application to condition estimators,
  507. //       SIAM J. Numer. Anal., 17 (1980), 403-409.
  508.  
  509. //    This file is a translation of qmult.m from version 2.0 of
  510. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  511. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  512.  
  513. //-------------------------------------------------------------------//
  514.  
  515. qmult = function ( A )
  516. {
  517.   local (A)
  518.  
  519.   n = A.nr; m = A.nc;
  520.  
  521.   //  Handle scalar A.
  522.   if (max(n,m) == 1)
  523.   {
  524.     n = A;
  525.     A = eye(n,n);
  526.   }
  527.  
  528.   d = zeros(n,n);
  529.  
  530.   for (k in n-1:1:-1)
  531.   {
  532.     // Generate random Householder transformation.
  533.     rand("normal", 0, 1);
  534.     x = rand(n-k+1,1);
  535.     s = norm(x, "2");
  536.     sgn = sign(x[1]) + (x[1]==0);    // Modification for sign(1)=1.
  537.     s = sgn*s;
  538.     d[k] = -sgn;
  539.     x[1] = x[1] + s;
  540.     beta = s*x[1];
  541.  
  542.     // Apply the transformation to A.
  543.     y = x'*A[k:n;];
  544.     A[k:n;] = A[k:n;] - x*(y/beta);
  545.   }
  546.  
  547.   // Tidy up signs.
  548.   for (i in 1:n-1)
  549.   {
  550.     A[i;] = d[i]*A[i;];
  551.   }
  552.  
  553.   A[n;] = A[n;]*sign(rand());
  554.   B = A;
  555.  
  556.   rand("default");
  557.   return B;
  558. };
  559.  
  560. //-------------------------------------------------------------------//
  561.  
  562. // Synopsis:    Create a Householder matrix.
  563.  
  564. // Syntax:    house ( X )
  565.  
  566. //    If HOUSE(x), which returns a list containing elements `v' and
  567. //    `beta',  then H = EYE - beta*v*v' is a Householder matrix such
  568. //    that Hx = -sign(x(1))*norm(x)*e_1. 
  569. //    NB: If x = 0 then v = 0, beta = 1 is returned.
  570. //            x can be real or complex.
  571. //            sign(x) := exp(i*arg(x)) ( = x./abs(x) when x ~= 0).
  572.  
  573. //    Theory: (textbook references Golub & Van Loan 1989, 38-43;
  574. //             Stewart 1973, 231-234, 262; Wilkinson 1965, 48-50).
  575. //      Hx = y: (I - beta*v*v')x = -s*e_1.
  576. //      Must have |s| = norm(x), v = x+s*e_1, and
  577. //      x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)).
  578. //      So take s = sign(x(1))*norm(x) (which avoids cancellation).
  579. //      v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2
  580. //          = 2*norm(x)*(norm(x) + |x(1)|).
  581.  
  582. //    References:
  583. //        G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  584. //           Johns Hopkins University Press, Baltimore, Maryland, 1989.
  585. //        G.W. Stewart, Introduction to Matrix Computations, Academic Press,
  586. //           New York, 1973,
  587. //        J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
  588. //           Press, 1965.
  589.  
  590. //    This file is a translation of house.m from version 2.0 of
  591. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  592. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  593.  
  594. //-------------------------------------------------------------------//
  595.  
  596. house = function ( x )
  597. {
  598.   local (x)
  599.  
  600.   m = x.nr; n = x.nc;
  601.   if (n > 1) { error ("Argument must be a column vector."); }
  602.  
  603.   s = norm(x,"2") * (sign(x[1]) + (x[1]==0)); // Modification for sign(0)=1.
  604.   v = x;
  605.   if (s == 0)    // Quit if x is the zero vector.
  606.   {
  607.     beta = 1; 
  608.     return << beta = beta ; v = v >>;
  609.   }
  610.  
  611.   v[1] = v[1] + s;
  612.   beta = 1/(s'*v[1]);        // NB the conjugated s.
  613.  
  614.   // beta = 1/(abs(s)*(abs(s)+abs(x(1)) would guarantee beta real.
  615.   // But beta as above can be non-real (due to rounding) only 
  616.   // when x is complex.
  617.  
  618.   return << beta = beta ; v = v >>;
  619. };
  620.  
  621. //-------------------------------------------------------------------//
  622.  
  623. //  Syntax:    tril ( A )
  624. //        tril ( A , K )
  625.  
  626. //  Description:
  627.  
  628. //  tril(x) returns the lower triangular part of A.
  629.  
  630. //  tril(A,K) returns the elements on and below the K-th diagonal of
  631. //  A.
  632.  
  633. //  K = 0: main diagonal
  634. //  K > 0: above the main diag.
  635. //  K < 0: below the main diag.
  636.  
  637. //  See Also: triu
  638. //-------------------------------------------------------------------//
  639.  
  640. tril = function(x, k) 
  641. {
  642.   local(x, k)
  643.  
  644.   if (!exist (k)) { k = 0; }
  645.   nr = x.nr; nc = x.nc;
  646.   if(k > 0) 
  647.   { 
  648.     if (k > (nc - 1)) { error ("tril: invalid value for k"); }
  649.   else
  650.     if (abs (k) > (nr - 1)) { error ("tril: invalid value for k"); }
  651.   }
  652.  
  653.   y = zeros(nr, nc);
  654.  
  655.   for(i in max( [1,1-k] ):nr) {
  656.     j = 1:min( [nc, i+k] );
  657.     y[i;j] = x[i;j];
  658.   }
  659.  
  660.   return y;
  661. };
  662.  
  663. //-------------------------------------------------------------------//
  664.  
  665. //  Syntax:    triu ( A )
  666. //        triu ( A , K )
  667.  
  668. //  Description:
  669.  
  670. //  triu(x) returns the upper triangular part of A.
  671.  
  672. //  tril(x; k) returns the elements on and above the k-th diagonal of
  673. //  A. 
  674.  
  675. //  K = 0: main diagonal
  676. //  K > 0: above the main diag.
  677. //  K < 0: below the main diag.
  678.  
  679. //  See Also: tril
  680. //-------------------------------------------------------------------//
  681.  
  682. triu = function(x, k) 
  683. {
  684.   local(x, k)
  685.  
  686.   if (!exist (k)) { k = 0; }
  687.   nr = x.nr; nc = x.nc;
  688.  
  689.   if(k > 0) 
  690.   { 
  691.     if (k > (nc - 1)) { error ("triu: invalid value for k"); }
  692.   else
  693.     if (abs (k) > (nr - 1)) { error ("triu: invalid value for k"); }
  694.   }
  695.  
  696.   y = zeros(nr, nc);
  697.  
  698.   for(j in max( [1,1+k] ):nc) {
  699.     i = 1:min( [nr, j-k] );
  700.     y[i;j] = x[i;j];
  701.   }
  702.  
  703.   return y;
  704. };
  705.  
  706. //
  707. //-------------------- Test Relational Expressions -------------------
  708. //
  709. printf("\tstart scalar tests...\n");
  710. printf("\tstart relational tests...\n");
  711.  
  712. //    SCALAR CONSTANTS (REAL)
  713. if( !(1<2) ) { error(); }
  714. if( !(1<=2) ) { error(); }
  715. if( 1>2 ) { error(); }
  716. if( 1>=2 ) { error(); }
  717. if( 1==2 ) { error(); }
  718. if( !(1!=2) ) { error(); }
  719. if( !1 ) { error(); }
  720. if( !!!1) { error(); }
  721.  
  722. if( !([1]<[2]) ) { error(); }
  723. if( !([1]<=[2]) ) { error(); }
  724. if( [1]>[2] ) { error(); }
  725. if( [1]>=[2] ) { error(); }
  726. if( [1]==[2] ) { error(); }
  727. if( !([1]!=[2]) ) { error(); }
  728. if( ![1] ) { error(); }
  729. if( !!![1]) { error(); }
  730.  
  731. //    SCALAR CONSTANTS (COMPLEX)
  732. if( !(1+2i<2+3i) ) { error(); }
  733. if( !(1+2i<=2+3i) ) { error(); }
  734. if( 1+2i>2+3i ) { error(); }
  735. if( 1+2i>=2+3i ) { error(); }
  736. if( 1+2i==2+3i ) { error(); }
  737. if( !(1+2i!=2+3i) ) { error(); }
  738. if( !1+2i ) { error(); }
  739. if( !!!1+2i) { error(); }
  740.  
  741. if( !([1+2i]<[2+3i]) ) { error(); }
  742. if( !([1+2i]<=[2+3i]) ) { error(); }
  743. if( [1+2i]>[2+3i] ) { error(); }
  744. if( [1+2i]>=[2+3i] ) { error(); }
  745. if( [1+2i]==[2+3i] ) { error(); }
  746. if( !([1+2i]!=[2+3i]) ) { error(); }
  747. if( ![1+2i] ) { error(); }
  748. if( !!![1+2i]) { error(); }
  749.  
  750. //    SCALAR ENTITIES (REAL)
  751. a=1;b=2;
  752. if( !(a<b) ) { error(); }
  753. if( !(a<=b) ) { error(); }
  754. if( a>b ) { error(); }
  755. if( a>=b ) { error(); }
  756. if( a==b ) { error(); }
  757. if( !(a!=b) ) { error(); }
  758. if( !a ) { error(); }
  759. if( !!!a) { error(); }
  760.  
  761. if( !([a]<[b]) ) { error(); }
  762. if( !([a]<=[b]) ) { error(); }
  763. if( [a]>[b] ) { error(); }
  764. if( [a]>=[b] ) { error(); }
  765. if( [a]==[b] ) { error(); }
  766. if( !([a]!=[b]) ) { error(); }
  767. if( ![a] ) { error(); }
  768. if( !!![a]) { error(); }
  769.  
  770. //    SCALAR ENTITIES (COMPLEX)
  771. a=1+2i;b=2+3i;
  772. if( !(a<b) ) { error(); }
  773. if( !(a<=b) ) { error(); }
  774. if( a>b ) { error(); }
  775. if( a>=b ) { error(); }
  776. if( a==b ) { error(); }
  777. if( !(a!=b) ) { error(); }
  778. if( !a ) { error(); }
  779. if( !!!a) { error(); }
  780.  
  781. if( !([a]<[b]) ) { error(); }
  782. if( !([a]<=[b]) ) { error(); }
  783. if( [a]>[b] ) { error(); }
  784. if( [a]>=[b] ) { error(); }
  785. if( [a]==[b] ) { error(); }
  786. if( !([a]!=[b]) ) { error(); }
  787. if( ![a] ) { error(); }
  788. if( !!![a]) { error(); }
  789.  
  790. x = rand(4,4);
  791.  
  792. if (! all (all (-2 <  x))) { "-2 < x" error (); }
  793. if (! all (all (-2 <= x))) { "-2 <= x" error (); }
  794. if (! all (all ( 2 >  x))) { "2 > x" error (); }
  795. if (! all (all ( 2 >= x))) { "2 >= x" error (); }
  796.  
  797. if (! all (all (rand(4,4) >  -2))) { error (); }
  798. if (! all (all (rand(4,4) >= -2))) { error (); }
  799. if (! all (all (rand(4,4) <   2))) { error (); }
  800. if (! all (all (rand(4,4) <=  2))) { error (); }
  801.  
  802. if (! all (all (rand (4,4) >  -rand (4,4)))) { error (); }
  803. if (! all (all (rand (4,4) >= -rand (4,4)))) { error (); }
  804. if (! all (all (-rand (4,4) <  rand (4,4)))) { error (); }
  805. if (! all (all (-rand (4,4) <= rand (4,4)))) { error (); }
  806.  
  807. //------------------------- Test REAL SCALARS ------------------------
  808. //
  809. //    CONSTANTS
  810. //      Addition
  811. if(1+2 != 3) { error(); }
  812. //      Subtraction
  813. if(1-2 != -1) { error(); }
  814. //      Multiply
  815. if(1*2 != 2) { error(); }
  816. //      Divide
  817. if(1/2 != 0.5) { error(); }
  818. //      Power
  819. if(2^2 != 4) { error(); }
  820. if(4^0 != 1) { error(); }
  821. //      Unary Minus
  822. if(2-3 != -1) { error(); }
  823. //
  824. //    ENTITIES
  825. //
  826. a = 1; b = 2; c = 3; d = 0.5;
  827. //      Addition
  828. if(a+b != c) { error(); }
  829. //      Subtraction
  830. if(a-b != -a) { error(); }
  831. //      Multiply
  832. if(a*b != b) { error(); }
  833. //      Divide
  834. if(a/b != d) { error(); }
  835. //      Power
  836. if(b^b != b*b) { error(); }
  837. if(b^0 != 1) { error (); }
  838. //      Unary Minus
  839. if(b-c != -a) { error(); }
  840. //
  841. //  ENTITIES & CONSTANTS
  842. //
  843. if(a+2 != c) { error(); }
  844. //      Subtraction
  845. if(1-b != -a) { error(); }
  846. //      Multiply
  847. if(a*2 != 2) { error(); }
  848. //      Divide
  849. if(1/b != d) { error(); }
  850. //      Power
  851. if(2^b != b*b) { error(); }
  852. //      Unary Minus
  853. if(b-3 != -a) { error(); }
  854. //
  855. //------------------------Test COMPLEX SCALARS -------------------------
  856. //
  857. //    CONSTANTS
  858. if(sqrt(-1) != 1i) { error(); }
  859. //      Addition
  860. if((1+2i)+(2+3i) != (3+5i)) { error(); }
  861. //      Subtraction
  862. if((1+2i)-(3+4i) != (-2-2i)) { error(); }
  863. //      Multiply
  864. if((1+2i)*(3+4i) != (-5+10i)) { error(); }
  865. //      Divide
  866. if((1+2i)/(3-4i) != (-.2+.4i)) { error(); }
  867. //      Power
  868. //      Precision problems prevent us from testing these. Have to
  869. //      be checked by hand.
  870. //  (1+2i)^2 = -3 + 4i
  871. //  (1+2i)^.5 = 1.272 + 7.862e-1i
  872. //  if((1+2i)^2 != (-3+4i)) { error(); }
  873. //  if((1+2i)^10 != (237+3116i)) { error(); }
  874. //      Unary Minus
  875. if(-(1+2i) != -1-2i) { error(); }
  876. //
  877. //    ENTITIES
  878. //
  879. a = 1+2i; b = 3+4i; c = 4+6i;
  880. if(a+b != c) { error(); }
  881. //      Subtraction
  882. if(a-b != -2-2i) { error(); }
  883. //      Multiply
  884. if(a*b != -5+10i) { error(); }
  885. //      Divide
  886. //if(a/(3-4i) != -.2+.4i) { error(); }
  887. //      Power
  888. //  if(b^b != b*b) { error(); }
  889. //      Unary Minus
  890. if(-a != -1-2i) { error(); }
  891. //
  892. //    ENTITIES & CONSTANTS
  893. //
  894. if(a+(3+4i) != c) { error(); }
  895. //      Subtraction
  896. if((1+2i)-b != -2-2i) { error(); }
  897. //      Multiply
  898. if(a*(3+4i) != -5+10i) { error(); }
  899. //      Divide
  900. //if(a/(3-4i) != -.2+.4i) { error(); }
  901. //      Power
  902. //if(b^b != b*b) { error(); }
  903. //      Unary Minus
  904. if(-(1+2i) != -1-2i) { error(); }
  905. //
  906. // String - Numerical Equalities
  907. //
  908.  
  909. if ((1 == "1")) { error(); }
  910. if (([1] == "1")) { error(); }
  911. if (("1" == 1)) { error(); }
  912. if (("1" == [1])) { error(); }
  913.  
  914. if (rand(3,3) == "str") { error(); }
  915.  
  916. if ("str" == rand(3,3)) { error(); }
  917. if (!any (any ((["1", "2"; "3", "4"] == "1") == [1,0;0,0]))) { error(); }
  918.  
  919. //
  920. //----------------------- Test REAL MATRICES ---------------------------
  921. //
  922.  
  923. printf("\tstart matrix tests...\n");
  924. printf("\t\treal-matrices\n");
  925.  
  926. //  Read in test matrices
  927. //
  928.  
  929. read("test_input");
  930.  
  931. //
  932. //  Matrix construction
  933. //
  934.  
  935. if(any([1;2;3] != [1,2,3]')) {
  936.   error();
  937. }
  938.  
  939. if(any (any (m0 != zeros(2,2)))) { error(); }
  940. if(any (any (m1 != 1+zeros(2,2)))) { error(); }
  941. if(any (any (m2 != [1,2;3,4]))) { error(); }
  942. if(any (any (m3 != [1+2i,2+3i;3+4i,5+6i]))) { error(); }
  943. if(any (any ([1,2;3+0i,4+0i] != m2))) { error(); }
  944. if(any (any ([m2] != m2))) { error(); }
  945.  
  946. //
  947. //  Matrix indexing
  948. //
  949.  
  950. p = pascal(6);
  951. if (!all (all (p[ [1:3] ; ] == p[ [1:3]' ; ]))) { error (); }
  952. if (!all (all (p[ ; [1:3] ] == p[ ; [1:3]' ]))) { error (); }
  953. if (!all (all (p[ [2:6] ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
  954. if (!all (all (p[ [2:6]' ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
  955. if (!all (all (p[ [6:12]' ] == p[ [6:12] ]))) { error (); }
  956.  
  957. //
  958. //  Sub-Matrix promotion
  959. //
  960.  
  961. if(m2[2;2] != 4) { error(); }
  962. if(any(m2[2;] != [3,4])) { error(); }
  963. if(any(m2[;2] != [2,4]')) { error(); }
  964. i=2;j=1;
  965. if(m2[i;j] != 3) { error(); }
  966. i=1;j=2;
  967. if(m2[i;j] != 2) { error(); }
  968. m = [1,2,3;4,5,6;7,8,9];
  969.  
  970. if(any(m[1;1,2] != [1,2])) 
  971. {
  972.   error();
  973. }
  974.  
  975. if(any(m[1,2;1] != [1;4])) 
  976. {
  977.   error();
  978. }
  979.  
  980. if(any (any (m[1,2;1,2] != [1,2;4,5]))) 
  981. {
  982.   error();
  983. }
  984.  
  985. //
  986.  
  987. if(m3[2;2] != (5+6i)) { error(); }
  988. if(any(m3[2;] != [3+4i,5+6i])) { error(); }
  989. if(any(m3[;2] != conj([2+3i,5+6i]'))) { error(); }
  990.  
  991. //
  992. //  Automatic creation, extension
  993. //
  994.  
  995. if(any (any ((mm[3;3]=10) != [0,0,0;0,0,0;0,0,10]))) { error(); }
  996. a=[1,2,3;4,5,6];
  997. if(any (any ((a[3;1]=10) != [1,2,3;4,5,6;10,0,0]))) { error(); }
  998. a=[1,2;3,4];
  999. if(any (any ((a[3,4;3,4]=[5,6;7,8]) != [1,2,0,0;3,4,0,0;0,0,5,6;0,0,7,8]))) 
  1000. {
  1001.   error();
  1002. }
  1003.  
  1004. mmm[2;] = a[1;];
  1005.  
  1006. //
  1007. //  Matrix binary operations
  1008. //
  1009.  
  1010. a = m2; b = [5,6;7,8];
  1011. if(any (any (a+a != [2,4;6,8]))) { error(); }
  1012. if(any (any (a-a != zeros(2,2)))) { error(); }
  1013. if(any (any (2+a != [3,4;5,6]))) { error(); }
  1014. if(any (any (2-a != [1,0;-1,-2]))) { error(); }
  1015. if(any (any (a-2 != [-1,0;1,2]))) { error(); }
  1016. if(any (any (2*a != [2,4;6,8]))) { error(); }
  1017. if(any (any ((a./2 != [0.5,1;1.5,2])))) { error(); }
  1018. if(any (any (a*a != [7,10;15,22]))) { error(); }
  1019. if(any (any (a*a*a != [37,54;81,118]))) { error(); }
  1020. if(any (any (a .* a != [1,4;9,16]))) { error(); }
  1021. if(any (any (a./a != [1,1;1,1]))) { error(); }
  1022. if(any (any (a' != [1,3;2,4]))) { error(); }
  1023.  
  1024. if(any(any(rand(3,3)^0 != eye(3,3)))) { error(); }
  1025. if(any(any(rand(3,3).^0 != ones(3,3)))) { error(); }
  1026. if(any(any(rand(1,3).^0 != ones(1,3)))) { error(); }
  1027. if(any(any(rand(3,1).^0 != ones(3,1)))) { error(); }
  1028. if(any(any(1.^zeros(3,1) != ones(3,1)))) { error(); }
  1029. if(any(any(1.^zeros(1,3) != ones(1,3)))) { error(); }
  1030.  
  1031. if(any ([1;2;3]+[4;5;6] != [5;7;9])) 
  1032. {
  1033.   error();
  1034. }
  1035.  
  1036. if(any ([1;2;3]-[4;5;6] != [-3;-3;-3])) 
  1037. {
  1038.   error();
  1039. }
  1040.  
  1041. if(any ([2;2;2] ./ [1;1;1] != [2;2;2])) 
  1042. {
  1043.   error();
  1044. }
  1045.  
  1046. if(any ([1;2;3] .* [4;5;6] != [4;10;18])) 
  1047. {
  1048.   error();
  1049. }
  1050.  
  1051. if (type (1^.33333) != "real") { error (); }
  1052. if (type (1^(1/3)) != "real") { error (); }
  1053. if (type ([1]^.33333) != "real") { error (); }
  1054. if (type (1^[.33333]) != "real") { error (); }
  1055. if (type ([1]^[.33333]) != "real") { error (); }
  1056.  
  1057. //
  1058. // Test row-wise matrix addition
  1059. //
  1060.  
  1061. a = [1,2,3]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1062. if (!all (all (a .+ b == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  1063. if (!all (all (b .+ a == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  1064.  
  1065. a = [1,2,3] + [1,2,3]*1i; 
  1066. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  1067. c = [2,4,6;5,7,9;8,10,12;11,13,15] + [2,4,6;5,7,9;8,10,12;11,13,15]*1i;
  1068. if (!all (all (a .+ b == c))) { error (); }
  1069. if (!all (all (b .+ a == c))) { error (); }
  1070.  
  1071. printf("\t\tpassed matrix row-wise add test...\n");
  1072.  
  1073. //
  1074. // Test row-wise matrix subtraction
  1075. //
  1076.  
  1077. a = [1,1,1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1078. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1079. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1080.  
  1081. a = [1,1,1] + [1,1,1]*1i;
  1082. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  1083. c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
  1084. if (!all (all (a .- b == -c))) { error (); }
  1085. if (!all (all (b .- a ==  c))) { error (); }
  1086.  
  1087. printf("\t\tpassed matrix row-wise subtraction test...\n");
  1088.  
  1089. //
  1090. // Test col-wise matrix addition
  1091. //
  1092.  
  1093. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1094. if (!all (all (a .+ b == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  1095. if (!all (all (b .+ a == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  1096.  
  1097. a = [1;1;1;1] + [1;1;1;1]*1i;
  1098. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  1099. c = [2,3,4;5,6,7;8,9,10;11,12,13] + [2,3,4;5,6,7;8,9,10;11,12,13]*1i;
  1100. if (!all (all (a .+ b == c))) { error (); }
  1101. if (!all (all (b .+ a == c))) { error (); }
  1102.  
  1103. printf("\t\tpassed matrix col-wise add test...\n");
  1104.  
  1105. //
  1106. // Test col-wise matrix subtraction
  1107. //
  1108.  
  1109. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1110. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1111. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1112.  
  1113. a = [1;1;1;1] + [1;1;1;1]*1i;
  1114. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  1115. c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
  1116. if (!all (all (a .- b == -c))) { error (); }
  1117. if (!all (all (b .- a ==  c))) { error (); }
  1118.  
  1119. printf("\t\tpassed matrix col-wise subtraction test...\n");
  1120.  
  1121. a = [1,2,3];
  1122. b = [1,2,3;4,5,6;7,8,9];
  1123. c = [1,4,9;4,10,18;7,16,27];
  1124.  
  1125. if (!all (all (a .* b == c))) { error (); }
  1126. if (!all (all (b .* a == c))) { error (); }
  1127.  
  1128. za = a + rand (size (a))*1j;
  1129. zb = b + rand (size (b))*1j;
  1130.  
  1131. if (!all (all (za .* zb == [za;za;za] .* zb))) { error (); }
  1132. if (!all (all (zb .* za == zb .* [za;za;za]))) { error (); }
  1133. if (!all (all (a .* zb == [a;a;a] .* zb))) { error (); }
  1134. if (!all (all (zb .* a == zb .* [a;a;a]))) { error (); }
  1135. if (!all (all (za .* b == [za;za;za] .* b))) { error (); }
  1136. if (!all (all (b .* za == b .* [za;za;za]))) { error (); }
  1137.  
  1138. printf("\t\tpassed matrix row-wise multiplication test...\n");
  1139.  
  1140. a = [1,2,3];
  1141. b = [1,2,3;4,6,6;7,8,9];
  1142. c = [1,1,1;4,3,2;7,4,3];
  1143.  
  1144. if (!all (all (b ./ a == c))) { error (); }
  1145. if (!all (all ([a;a;a] ./ b == a ./ b))) { error (); }
  1146. if (!all (all (b ./ [a;a;a] == b ./ a))) { error (); }
  1147.  
  1148. za = a + rand (size (a))*1j;
  1149. zb = b + rand (size (b))*1j;
  1150.  
  1151. if (!all (all ([za;za;za] ./ zb == za ./ zb))) { error (); }
  1152. if (!all (all (zb ./ [za;za;za] == zb ./ za))) { error (); }
  1153. if (!all (all ([a;a;a] ./ zb == a ./ zb))) { error (); }
  1154. if (!all (all (zb ./ [a;a;a] == zb ./ a))) { error (); }
  1155. if (!all (all ([za;za;za] ./ b == za ./ b))) { error (); }
  1156. if (!all (all (b ./ [za;za;za] == b ./ za))) { error (); }
  1157.  
  1158. printf("\t\tpassed matrix row-wise division test...\n");
  1159.  
  1160. a = [1;2;3];
  1161. b = [1,2,3;4,5,6;7,8,9];
  1162.  
  1163. if (!all (all (a .* b == [a,a,a] .* b))) { error (); }
  1164. if (!all (all (b .* a == b .* [a,a,a]))) { error (); }
  1165.  
  1166. za = a + rand (size (a))*1j;
  1167. zb = b + rand (size (b))*1j;
  1168.  
  1169. if (!all (all (za .* zb == [za,za,za] .* zb))) { error (); }
  1170. if (!all (all (zb .* za == zb .* [za,za,za]))) { error (); }
  1171. if (!all (all (za .* b == [za,za,za] .* b))) { error (); }
  1172. if (!all (all (b .* za == b .* [za,za,za]))) { error (); }
  1173. if (!all (all (a .* zb == [a,a,a] .* zb))) { error (); }
  1174. if (!all (all (zb .* a == zb .* [a,a,a]))) { error (); }
  1175.  
  1176. printf("\t\tpassed matrix column-wise multiplication test...\n");
  1177.  
  1178. a = [1;2;3];
  1179. b = [1,2,3;4,6,6;7,8,9];
  1180.  
  1181. if (!all (all ([a,a,a] ./ b == a ./ b))) { error (); }
  1182. if (!all (all (b ./ [a,a,a] == b ./ a))) { error (); }
  1183.  
  1184. za = a + rand (size(a))*1j;
  1185. zb = b + rand (size(b))*1j;
  1186.  
  1187. if (!all (all ([za,za,za] ./ zb == za ./ zb))) { error (); }
  1188. if (!all (all (zb ./ [za,za,za] == zb ./ za))) { error (); }
  1189. if (!all (all ([za,za,za] ./ b == za ./ b))) { error (); }
  1190. if (!all (all (b ./ [za,za,za] == b ./ za))) { error (); }
  1191. if (!all (all ([a,a,a] ./ zb == a ./ zb))) { error (); }
  1192. if (!all (all (zb ./ [a,a,a] == zb ./ a))) { error (); }
  1193.  
  1194. printf("\t\tpassed matrix column-wise division test...\n");
  1195.  
  1196.  
  1197. //
  1198. //--------------------- Test COMPLEX MATRICES --------------------------
  1199. //
  1200. //  Automatic creation, extension
  1201. //
  1202. printf("\t\tcomplex-matrices\n");
  1203. if(any (any ((mm[3;3]=10+10i) != [0,0,0;0,0,0;0,0,10+10i]))) { error(); }
  1204. a=[1,2,3;4,5,6];
  1205. if(any (any ((a[3;1]=10+10i) != [1,2,3;4,5,6;10+10i,0,0]))) { error(); }
  1206. //
  1207. a = m3;
  1208. if(any (any (a+a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  1209. if(any (any (a-a != zeros(2,2)))) { error(); }
  1210. if(any (any (2+a != [3+2i,4+3i;5+4i,7+6i]))) { "2+a" error(); }
  1211. if(any (any (2-a != [1-2i,0-3i;-1-4i,-3-6i]))) { error(); }
  1212. if(any (any (a-2 != [-1+2i,0+3i;1+4i,3+6i]))) { error(); }
  1213. if(any (any (2*a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  1214. if(any (any (a./2 != [.5+1i,1+1.5i;1.5+2i,2.5+3i]))) { error(); }
  1215. if(any (any (a*a != [-9+21i,-12+34i;-14+48i,-17+77i]))) { error(); }
  1216. if(any (any (a*a*a != [-223+57i,-345+113i;-469+183i,-719+337i]))) { error(); }
  1217. if(any (any (a .* a != [-3+4i,-5+12i;-7+24i,-11+60i]))) { error(); }
  1218. //
  1219. // The following test may not work on some machines
  1220. //
  1221. if(any (any (a./a != [1,1;1,1]))) { 
  1222.   printf("\t\t***complex division inaccuracy, check manually***\n");
  1223. }
  1224.  
  1225. if(any (any (a' != conj([1+2i,3+4i;2+3i,5+6i])))) { error(); }
  1226. //  
  1227. //--------------------- Test NULL MATRICES -------------------------
  1228. //
  1229. printf("\t\tnull-matrices\n");
  1230. // Create a NULL matrix
  1231. mnull = [];
  1232. // Test it with SCALARS
  1233. if( any([1,mnull] != 1)) {
  1234.   error();
  1235. }
  1236. if( any([mnull,1] != 1)) {
  1237.   error();
  1238. }
  1239. // Test with MATRIX construction
  1240. m = [1,2;3,4;5,6];
  1241. if( any([mnull;1] != [1])) {
  1242.   error();
  1243. }
  1244. if( any([1;mnull] != [1])) {
  1245.   error();
  1246. }
  1247. if( any([mnull;1,2,3] != [1,2,3])) {
  1248.   error();
  1249. }
  1250. if( any([1,2,3;mnull] != [1,2,3])) {
  1251.   error();
  1252. }
  1253. if(any (any ([mnull,m] != m))) {
  1254.   error();
  1255. }
  1256. if(any (any ([m,mnull] != m))) {
  1257.   error();
  1258. }
  1259. if(any (any ([mnull;m] != m))) {
  1260.   error();
  1261. }
  1262. if(any (any ([m;mnull] != m))) {
  1263.   error();
  1264.  
  1265. mnull = matrix();
  1266. mnull[1] = [1];
  1267. }
  1268.  
  1269. //--------------------- Test Matrix Multiply --------------------------
  1270.  
  1271. i = sqrt(-1);
  1272. a = [1,2,3;4,5,6;7,8,9];
  1273. b = [4,5,6;7,8,9;10,11,12];
  1274. c = [ 48,  54,  60 ;
  1275.      111, 126, 141 ;
  1276.      174, 198, 222 ];
  1277.  
  1278. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  1279.  
  1280. az = a + b*i;
  1281. bz = b + a*i;
  1282.  
  1283. cz = [-18+141*i , -27+162*i , -36+183*i ;
  1284.         9+240*i ,   0+279*i ,  -9+318*i ;
  1285.        36+339*i ,  27+396*i ,  18+453*i ];
  1286.  
  1287. czz = [ 48+30*i ,  54+36*i  ,  60+42*i ;
  1288.        111+66*i , 126+81*i  , 141+96*i ;
  1289.        174+102*i, 198+126*i , 222+150*i ];
  1290.  
  1291. czzz = [ 48+111*i ,  54+126*i ,  60+141*i ;
  1292.         111+174*i , 126+198*i , 141+222*i ;
  1293.         174+237*i , 198+270*i , 222+303*i ];
  1294.  
  1295. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  1296. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  1297. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  1298.  
  1299. a = [a,a];
  1300. b = [b;b];
  1301. c = [  96 , 108 , 120 ;
  1302.       222 , 252 , 282 ;
  1303.       348 , 396 , 444 ];
  1304.  
  1305. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  1306.  
  1307. az = [az,az];
  1308. bz = [bz;bz];
  1309.  
  1310. cz = [  -36+282*i ,  -54+324*i ,  -72+366*i ;
  1311.          18+480*i ,    0+558*i ,  -18+636*i ;
  1312.          72+678*i ,   54+792*i ,   36+906*i ];
  1313.  
  1314. czz = [  96+60*i  , 108+72*i  , 120+84*i  ;
  1315.         222+132*i , 252+162*i , 282+192*i ;
  1316.         348+204*i , 396+252*i , 444+300*i ];
  1317.  
  1318. czzz = [  96+222*i , 108+252*i , 120+282*i ;
  1319.          222+348*i , 252+396*i , 282+444*i ;
  1320.          348+474*i , 396+540*i , 444+606*i ];
  1321.  
  1322. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  1323. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  1324. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  1325.  
  1326. printf("\t\tpassed matrix multiply test...\n");
  1327.  
  1328. //--------------------- Test STRING MATRICES --------------------------
  1329. //
  1330. printf("\t\tstring-matrices\n");
  1331. sm = ["s1","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"];
  1332. if(sm[1] != "s1") { error(); }
  1333. if( sm[1;3] != "sm3" ) { error(); }
  1334. if(any(sm[2,3;3] != ["xxx3";"yyy3"]) ) { error(); }
  1335. if(any (any ((sm[1;1]="xx")!=["xx","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"])))
  1336. {
  1337.   error();
  1338. }
  1339. if( ((strm[1] = "strm")[1]) != "strm" ) { error(); }
  1340.  
  1341. //  Test string-matrix add functionality
  1342.  
  1343. sm = sm[1,2;1,2];
  1344. if (any (any (("1_"+sm+"_2") != ["1_xx_2","1_sm2_2";"1_x1_2","1_x2_2"]))) {error();}
  1345.  
  1346. if ("c"+["1"] != "c1") { error (); }
  1347. if (["c"]+"1" != "c1") { error (); }
  1348.  
  1349. // An append test...
  1350. bs = "b";
  1351. bm = "b";
  1352.  
  1353. xsm = ["a", bm, bs ];
  1354. if (!all(xsm == ["a", "b", "b"])) { error (); }
  1355.  
  1356. printf("\tpassed matrix test...\n");
  1357.  
  1358. //
  1359. //---------------------------- List Tests --------------------------
  1360. //
  1361. //  List creation
  1362. listest = << << 11; 12 >>; << 21; 22>> >>;
  1363. if( listest.[1].[2] != 12 ) { error(); }
  1364. if(any(<<a=10;b=1:4;c=[1,2,3;4,5,6;7,8,9]>>.b != [1,2,3,4])) { error(); }
  1365. mlist.[0] = m;
  1366. if(any(any(mlist.[0] != m))) { error(); }
  1367.  
  1368. // Test list functions...
  1369.  
  1370. listtest.fun = function ( a ) { return sin(a); };
  1371. if (!(listtest.fun(0.5) == sin(0.5))) { error() ; }
  1372.  
  1373. // Test open-list assignment...
  1374. </v;d/> = eig(rand(3,3));
  1375.  
  1376. printf("\tpassed list test...\n");
  1377.  
  1378. //
  1379. // Reset random number generator seed...
  1380. //
  1381.  
  1382. rand("default");
  1383.  
  1384. //
  1385. //-------------------------- Test printf () --------------------------
  1386. //
  1387.  
  1388. sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", 5,3,2, 8,7,3, "string", 3, 4, 1234e-2);
  1389. if (!(tmp == "  002  0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
  1390.  
  1391. sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", [5],[3],[2], [8],[7],[3], ...
  1392.          ["string"], [3], [4], [1234e-2]);
  1393. if (!(tmp == "  002  0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
  1394.  
  1395. sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", [5,1],[3,2],[2,2], [8],[7],[3], ...
  1396.          ["string"], [3], [4], [1234e-2,4]);
  1397. if (!(tmp == "  002  0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
  1398.  
  1399. sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", [5+2i,1],[3,2+4i],[2,2], [8],[7],[3], ...
  1400.          ["string"], [3+2i], [4], [1234e-2+12j,4]);
  1401. if (!(tmp == "  002  0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
  1402.  
  1403. printf("\tpassed sprintf test...\n");
  1404.  
  1405. //
  1406. //------------------------- Test strtod()  ----------------------------
  1407. //
  1408.  
  1409. if (123.456 != strtod ("123.456")) { error (); }
  1410. if (!all (all ([1,2;3,4] == strtod (["1","2";"3","4"]))))
  1411.   { error (); }
  1412. printf("\tpassed strtod test...\n");
  1413.  
  1414. //
  1415. //------------------------- Test getline()  ---------------------------
  1416. //
  1417. //
  1418.  
  1419. close( "test_getl" );
  1420.  
  1421. x = getline( "test_getl" );
  1422. if ( x.[1] !=  123.456 ) { error(); }
  1423. if ( x.[2] != -123.456 ) { error(); }
  1424. if ( x.[3] !=  123.456 ) { error(); }
  1425.  
  1426. x = getline( "test_getl" );
  1427. if ( x.[1] !=  .123 ) { error(); }
  1428. if ( x.[2] != -.123 ) { error(); }
  1429. if ( x.[3] !=  .123 ) { error(); }
  1430.  
  1431. x = getline( "test_getl" );
  1432. if ( x.[1] !=  123 ) { error(); }
  1433. if ( x.[2] != -123 ) { error(); }
  1434. if ( x.[3] !=  123 ) { error(); }
  1435.  
  1436. x = getline( "test_getl" );
  1437. if ( x.[1] !=  1e6 ) { error(); }
  1438. if ( x.[2] != -1e6 ) { error(); }
  1439. if ( x.[3] !=  1e6 ) { error(); }
  1440.  
  1441. x = getline( "test_getl" );
  1442. if ( x.[1] !=  1e5 ) { error(); }
  1443. if ( x.[2] != -1e5 ) { error(); }
  1444. if ( x.[3] !=  1e5 ) { error(); }
  1445.  
  1446. x = getline( "test_getl" );
  1447. if ( x.[1] !=  123.456e3 ) { error(); }
  1448. if ( x.[2] != -123.456e3 ) { error(); }
  1449. if ( x.[3] !=  123.456e3 ) { error(); }
  1450.  
  1451. x = getline( "test_getl" );
  1452. if ( x.[1] !=  123.456e3 ) { error(); }
  1453. if ( x.[2] != -123.456e3 ) { error(); }
  1454. if ( x.[3] !=  123.456e3 ) { error(); }
  1455.  
  1456. x = getline( "test_getl" );
  1457. if ( x.[1] !=  123.456e-3 ) { error(); }
  1458. if ( x.[2] != -123.456e-3 ) { error(); }
  1459. if ( x.[3] !=  123.456e-3 ) { error(); }
  1460.  
  1461. x = getline( "test_getl" );
  1462. if ( x.[1] !=  .123e3 ) { error(); }
  1463. if ( x.[2] != -.123e3 ) { error(); }
  1464. if ( x.[3] !=  .123e3 ) { error(); }
  1465.  
  1466. x = getline( "test_getl" );
  1467. if ( x.[1] !=  123e3 ) { error(); }
  1468. if ( x.[2] != -123e3 ) { error(); }
  1469. if ( x.[3] !=  123e3 ) { error(); }
  1470.  
  1471. x = getline( "test_getl" );
  1472. if ( x.[1] != "123abc" ) { error(); }
  1473. if ( x.[2] != "abc123" ) { error(); }
  1474. if ( x.[3] != "123.abc" ) { error(); }
  1475.  
  1476. x = getline( "test_getl" );
  1477. if ( x.[1] != "quoted string" ) { error(); }
  1478. if ( x.[2] != "q string with escapes \n \t \" " ) { error(); }
  1479.  
  1480. x = getline( "test_getl" );
  1481. if ( x.[1] != "quoted string" ) { error(); }
  1482. if ( x.[2] !=  1.23e3 ) { error(); }
  1483. if ( x.[3] !=  100 ) { error(); }
  1484. if ( x.[4] != "q string with escapes \n \t \" " ) { error(); }
  1485. if ( x.[5] !=  200 ) { error(); }
  1486.  
  1487. x = getline ("test_getl", 0);
  1488. if (type (x) != "string") { error (); }
  1489.  
  1490. // Also check out strsplt while we are here...
  1491. if (length (strsplt(x)) != 79) { error (); }
  1492. if (length (strsplt(x, 13)) != 6) { error (); }
  1493. if (length (strsplt(x,[13,12,14,13,15,11] )) != 6) { error (); }
  1494. if (length (strsplt(x, ".")) != 7) { error (); }
  1495.  
  1496. printf("\tpassed getline() test...\n");
  1497.  
  1498. //
  1499. //---------------------- Test readb()/writeb()  --------------------
  1500. //
  1501.  
  1502. a = rand (5,5);
  1503. z = rand(3,3) + rand(3,3)*1j;
  1504. strm = what()[1:5;1:5];
  1505. pi = 4*atan(1.0);
  1506. sc = 2*pi;
  1507. str = "this is a sample string\ttab";
  1508. l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1509.  
  1510. writeb ("jnk_rb", a, z, strm, sc, str, l);
  1511.  
  1512. # Set aside matrices for later tests
  1513.  
  1514. check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1515.  
  1516. clear (a, z, strm, sc, str, l);
  1517.  
  1518. close ("jnk_rb");
  1519.  
  1520. readb ("jnk_rb");
  1521.  
  1522. #
  1523. # Now do checks
  1524. #
  1525.  
  1526. if (!all (all (a == check.a))) { error (); }
  1527. if (!all (all (z == check.z))) { error (); }
  1528. if (!all (all (strm == check.strm))) { error (); }
  1529. if (sc != check.sc) { error (); }
  1530. if (str != check.str) { error (); }
  1531.  
  1532. if (length (l) != 5) { error (); }
  1533. if (!all (all (l.a == check.a))) { error (); }
  1534. if (!all (all (l.z == check.z))) { error (); }
  1535. if (!all (all (l.strm == check.strm))) { error (); }
  1536. if (l.sc != check.sc) { error (); }
  1537. if (l.str != check.str) { error (); }
  1538.  
  1539.  
  1540. printf("\tpassed binary I/O test...\n");
  1541.  
  1542. //
  1543. //---------------------- Test read()/write()  --------------------
  1544. //
  1545.  
  1546. a = rand (5,5);
  1547. z = rand(3,3) + rand(3,3)*1j;
  1548. strm = what()[1:5;1:5];
  1549. pi = 4*atan(1.0);
  1550. sc = 2*pi;
  1551. str = "this is a sample string\ttab";
  1552. l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1553.  
  1554. write ("jnk_ra", a, z, strm, sc, str, l);
  1555.  
  1556. # Set aside matrices for later tests
  1557.  
  1558. check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1559.  
  1560. clear (a, z, strm, sc, str, l);
  1561.  
  1562. close ("jnk_ra");
  1563.  
  1564. read ("jnk_ra");
  1565.  
  1566. #
  1567. # Now do checks
  1568. #
  1569.  
  1570. if (a.nr != 5 || a.nc != 5) { error (); }
  1571. if (z.nr != 3 || z.nc != 3) { error (); }
  1572. if (strm.nr != 5 || strm.nc != 5) { error (); }
  1573. if (str != check.str) { error (); }
  1574.  
  1575. if (length (l) != 5) { error (); }
  1576. if (l.a.nr != 5 || l.a.nc != 5) { error (); }
  1577. if (l.z.nr != 3 || l.z.nc != 3) { error (); }
  1578. if (l.strm.nr != 5 || l.strm.nc != 5) { error (); }
  1579. if (l.str != check.str) { error (); }
  1580.  
  1581. printf("\tpassed ascii I/O test...\n");
  1582.  
  1583. //
  1584. //-------------------------- Test eval () ------------------------------
  1585. //
  1586.  
  1587. if (1 + 2 != eval("1 + 2")) { error ("eval() error"); }
  1588. x = function (s, a) { return eval(s); };
  1589. str = "yy = 2 + x(\"2*a\",3)";
  1590. z = eval(str);
  1591. if (z != yy) { error ("eval() error"); }
  1592. printf("\tpassed eval test...\n");
  1593.  
  1594.  
  1595. //-------------------------- ------------------------------
  1596.  
  1597. printf ("\ttest builtin function correctness...\n");
  1598.  
  1599. //-------------------------- ------------------------------
  1600.  
  1601. //
  1602. //-------------------------- Test abs () ------------------------------
  1603. //
  1604.  
  1605. A = rand(5,5);
  1606. T = ( A == abs (A));
  1607. if (!all (all (A))) { error ("abs() incorrect"); }
  1608. printf("\tabs()");
  1609.  
  1610.  
  1611. //
  1612. //-------------------------- Test max () ------------------------------
  1613. //
  1614.  
  1615. A = [1,10,100;2,20,200;1,2,3];
  1616. B = A./2;
  1617. ZA = A + rand (3,3)*A*1i;
  1618. ZB = B + rand (3,3)*B*1i;
  1619. if (!all (max (A) == [2,20,200])) { error( "max() incorrect"); }
  1620. if (max (max(A)) != 200) { error ("max() incorrect"); }
  1621. if (any (any (max (A, B) != A))) { error (); }
  1622. if (any (any (max (B, A) != max (A, B)))) { error (); }
  1623. if (any (any (max (ZB, ZA) != max (ZA, ZB)))) { error (); }
  1624. if (any (any (max (B, ZA) != max (ZA, B)))) { error (); }
  1625. if (any (any (max (ZB, A) != max (A, ZB)))) { error (); }
  1626. printf("\tmax()");
  1627.  
  1628. //
  1629. //-------------------------- Test min () ------------------------------
  1630. //
  1631.  
  1632. if (!all (min (A) == [1,2,3])) { error( "min() incorrect"); }
  1633. if (min (min(A)) != 1) { error ("min() incorrect"); }
  1634. if (any (any (min (A, B) != B))) { error (); }
  1635. if (any (any (min (B, A) != min (A, B)))) { error (); }
  1636. if (any (any (min (ZB, ZA) != min (ZA, ZB)))) { error (); }
  1637. if (any (any (min (B, ZA) != min (ZA, B)))) { error (); }
  1638. if (any (any (min (ZB, A) != min (A, ZB)))) { error (); }
  1639. printf("\tmin()");
  1640.  
  1641. //
  1642. //-------------------------- Test maxi () -----------------------------
  1643. //
  1644.  
  1645. if (!all (maxi (A) == [2,2,2])) { error( "maxi() incorrect"); }
  1646. if (maxi (maxi(A)) != 1) { error ("maxi() incorrect"); }
  1647. printf("\tmaxi()");
  1648.  
  1649. //
  1650. //-------------------------- Test mini () -----------------------------
  1651. //
  1652.  
  1653. if (!all (mini (A) == [1,3,3])) { error( "mini() incorrect"); }
  1654. if (mini (mini(A)) != 1) { error ("mini() incorrect"); }
  1655. printf("\tmini()");
  1656.  
  1657. //
  1658. //-------------------------- Test floor () ----------------------------
  1659. //
  1660.  
  1661. if (floor (1.9999) != 1) { error ("floor() output incorrect"); }
  1662. if (!all (all (floor ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  1663.   error ("floor output incorrect");
  1664. }
  1665. printf("\tfloor()");
  1666.  
  1667. //
  1668. //-------------------------- Test ceil () 0----------------------------
  1669. //
  1670.  
  1671. if (ceil (1.9999) != 2) { error ("ceil() output incorrect"); }
  1672. if (!all (all (ceil ([1.99,1.99;2.99,2.99]) == [2,2;3,3]))) {
  1673.   error ("ceil output incorrect");
  1674. }
  1675. printf("\tceil()\n");
  1676.  
  1677. //
  1678. //-------------------------- Test round () ----------------------------
  1679. //
  1680.  
  1681. if (round (1.8) != 2) { error ("round() output incorrect"); }
  1682. if (round (1.4) != 1) { error ("round() output incorrect"); }
  1683. if (!all (all (round ([1.99,1.99;2.4,2.4]) == [2,2;2,2]))) {
  1684.   error ("round output incorrect");
  1685. }
  1686. printf("\tround()");
  1687.  
  1688. //
  1689. //-------------------------- Test sum () ------------------------------
  1690. //
  1691.  
  1692. S = [1:4; 4:7; 8:11];
  1693. if (sum (S[1;]) != 10) { error ("sum() incorrect"); }
  1694. if (sum (S[3;]) != 38) { error ("sum() incorrect"); }
  1695. if (!all (all (sum (S) == [13,16,19,22]))) { error ("sum() incorrect"); }
  1696. printf("\tsum()");
  1697.  
  1698.  
  1699. //
  1700. //-------------------------- Test cumsum () ------------------------------
  1701. //
  1702.  
  1703. S = [1:4; 4:7; 8:11];
  1704. if (any(cumsum (S[1;]) != [1,3,6,10])) { error ("cumsum() incorrect"); }
  1705. if (any(cumsum (S[;3]) != [3;9;19])) { error ("cumsum() incorrect"); }
  1706. if (!all (all (cumsum (S) == [1,2,3,4;5,7,9,11;13,16,19,22]))) 
  1707. { 
  1708.   error ("cumsum() incorrect"); 
  1709. }
  1710. printf("\tcumsum()");
  1711.  
  1712. //
  1713. //-------------------------- Test prod () ------------------------------
  1714. //
  1715.  
  1716. S = [1:4; 4:7; 8:11];
  1717. if (prod (S[1;]) != 24) { error ("prod() incorrect"); }
  1718. if (prod (S[;3]) != 180) { error ("prod() incorrect"); }
  1719. if (!all (all (prod (S) == [32,90,180,308]))) 
  1720. { 
  1721.   error ("prod() incorrect"); 
  1722. }
  1723. printf("\tprod()");
  1724.  
  1725. //
  1726. //-------------------------- Test cumprod () ------------------------------
  1727. //
  1728.  
  1729. S = [1:4; 4:7; 8:11];
  1730. if (any(cumprod (S[1;]) != [1,2,6,24])) { error ("cumprod() incorrect"); }
  1731. if (any(cumprod (S[;3]) != [3;18;180])) { error ("cumprod() incorrect"); }
  1732. if (!all (all (cumprod (S) == [1,2,3,4;4,10,18,28;32,90,180,308]))) 
  1733. { 
  1734.   error ("cumprod() incorrect"); 
  1735. }
  1736. printf("\tcumprod()\n");
  1737.  
  1738. //
  1739. //-------------------------- Test int () ------------------------------
  1740. //
  1741.  
  1742. if (int (1.9999) != 1) { error ("int() output incorrect"); }
  1743. if (!all (all (int ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  1744.   error ("int() output incorrect");
  1745. }
  1746. printf("\tint()");
  1747.  
  1748. //
  1749. //-------------------------- Test mod () ------------------------------
  1750. //
  1751.  
  1752. if (mod (1,1) != 0) { error ("mod() output incorrect"); }
  1753. if (mod (4,2) != 0) { error ("mod() output incorrect"); }
  1754. if (mod (3,2) != 1) { error ("mod() output incorrect"); }
  1755. if (mod (5,3) != 2) { error ("mod() output incorrect"); }
  1756. printf("\tmod()");
  1757.  
  1758. //
  1759. //-------------------------- Test find () ------------------------------
  1760. //
  1761.  
  1762. if (find ([0,1]) != 2) { error ("find() output incorrect"); }
  1763. if (find ([1,0]) != 1) { error ("find() output incorrect"); }
  1764. if (find ([0,1+1i]) != 2) { error ("find() output incorrect"); }
  1765. if (find ([1+0i,0]) != 1) { error ("find() output incorrect"); }
  1766. if (find ([0+1i,0]) != 1) { error ("find() output incorrect"); }
  1767.  
  1768. printf("\tfind()\n");
  1769.  
  1770.  
  1771. //-------------------------- ------------------------------
  1772.  
  1773. printf ("\ttest builtin function operation...\n");
  1774.  
  1775. //-------------------------- ------------------------------
  1776.  
  1777. //
  1778. // At least use most of the builtins....
  1779. //
  1780.  
  1781. A = rand(5,5);
  1782. Z = rand(5,5) + rand(5,5)*1j;
  1783. S = what()[3:6;2:4];
  1784.  
  1785. abs (A);
  1786. abs ([A]);
  1787. abs (3.14);
  1788.  
  1789. abs (Z);
  1790. abs ([Z]);
  1791. abs (3.14j);
  1792.  
  1793. printf ("\tabs");
  1794.  
  1795. acos (A);
  1796. acos ([A]);
  1797. acos (3.14);
  1798.  
  1799. acos (Z);
  1800. acos ([Z]);
  1801. acos (3.14j);
  1802.  
  1803. printf ("\tacos");
  1804.  
  1805. all (A);
  1806. all ([A]);
  1807. all (3.14);
  1808.  
  1809. all (Z);
  1810. all ([Z]);
  1811. all (3.14j);
  1812.  
  1813. printf ("\tall");
  1814.  
  1815. any (A);
  1816. any ([A]);
  1817. any (3.14);
  1818.  
  1819. any (Z);
  1820. any ([Z]);
  1821. any (3.14j);
  1822.  
  1823. printf ("\tany");
  1824.  
  1825. asin (A);
  1826. asin ([A]);
  1827. asin (3.14);
  1828.  
  1829. asin (Z);
  1830. asin ([Z]);
  1831. asin (3.14j);
  1832.  
  1833. printf ("\tasin");
  1834.  
  1835. atan (A);
  1836. atan ([A]);
  1837. atan (3.14);
  1838.  
  1839. atan (Z);
  1840. atan ([Z]);
  1841. atan (3.14j);
  1842.  
  1843. printf ("\tatan");
  1844.  
  1845. atan2 (A,A);
  1846. atan2 ([A],[A]);
  1847. atan2 (3.14,3.14);
  1848.  
  1849. #atan2 (Z,Z);
  1850. #atan2 ([Z],[Z]);
  1851. #atan2 (3.14j,3.14j);
  1852.  
  1853. printf ("\tatan2");
  1854. printf ("\n");
  1855.  
  1856. cd ("<RLaB$dir>");
  1857. printf ("\tcd");
  1858.  
  1859. ceil (A);
  1860. ceil ([A]);
  1861. ceil (3.14);
  1862.  
  1863. ceil (Z);
  1864. ceil ([Z]);
  1865. ceil (3.14j);
  1866.  
  1867. printf ("\tceil");
  1868.  
  1869. class (A);
  1870. class ([A]);
  1871. class (3.14);
  1872. class ("string");
  1873. class (S);
  1874. class ([S]);
  1875. class (<< rand(3,3); rand(4,4) >>);
  1876.  
  1877. class (Z);
  1878. class ([Z]);
  1879. class (3.14j);
  1880.  
  1881. printf ("\tclass");
  1882.  
  1883. clear (rand(3,3));
  1884. clear ([rand(3,3)]);
  1885.  
  1886. printf ("\tclear");
  1887.  
  1888. conj (A);
  1889. conj ([A]);
  1890. conj (3.14);
  1891.  
  1892. conj (Z);
  1893. conj ([Z]);
  1894. conj (3.14j);
  1895.  
  1896. printf ("\tconj");
  1897.  
  1898. cos (A);
  1899. cos ([A]);
  1900. cos (3.14);
  1901.  
  1902. cos (Z);
  1903. cos ([Z]);
  1904. cos (3.14j);
  1905.  
  1906. printf ("\tcos");
  1907.  
  1908. diag (A);
  1909. diag ([A]);
  1910. diag (3.14);
  1911.  
  1912. diag (Z);
  1913. diag ([Z]);
  1914. diag (3.14j);
  1915.  
  1916. printf ("\tdiag");
  1917. printf ("\n");
  1918.  
  1919. exist (A);
  1920. exist (Z);
  1921.  
  1922. printf ("\texist");
  1923.  
  1924. exp (A);
  1925. exp ([A]);
  1926. exp (3.14);
  1927.  
  1928. exp (Z);
  1929. exp ([Z]);
  1930. exp (3.14j);
  1931.  
  1932. printf ("\texp");
  1933.  
  1934. find (A);
  1935. find ([A]);
  1936. find (3.14);
  1937.  
  1938. find (Z);
  1939. find ([Z]);
  1940. find (3.14j);
  1941.  
  1942. printf ("\tfind");
  1943.  
  1944. floor (A);
  1945. floor ([A]);
  1946. floor (3.14);
  1947.  
  1948. floor (Z);
  1949. floor ([Z]);
  1950. floor (3.14j);
  1951.  
  1952. printf ("\tfloor");
  1953.  
  1954. format (1,1);
  1955. format ([2],[3]);
  1956. format ();
  1957. format ([3], 3);
  1958. format ();
  1959.  
  1960. printf ("\tformat");
  1961.  
  1962. imag (A);
  1963. imag ([A]);
  1964. imag (3.14);
  1965.  
  1966. imag (Z);
  1967. imag ([Z]);
  1968. imag (3.14j);
  1969.  
  1970. printf ("\timag");
  1971.  
  1972. inf ();
  1973. printf ("\tinf");
  1974. printf ("\n");
  1975.  
  1976. int (A);
  1977. int ([A]);
  1978. int (3.14);
  1979.  
  1980. int (Z);
  1981. int ([Z]);
  1982. int (3.14j);
  1983.  
  1984. printf ("\tint");
  1985.  
  1986. issymm (A);
  1987. issymm ([A]);
  1988. issymm (3.14);
  1989.  
  1990. issymm (Z);
  1991. issymm ([Z]);
  1992. issymm (3.14j);
  1993.  
  1994. printf ("\tissymm");
  1995.  
  1996. length (A);
  1997. length ([A]);
  1998. length (3.14);
  1999.  
  2000. length (Z);
  2001. length ([Z]);
  2002. length (3.14j);
  2003.  
  2004. length (3.14);
  2005. length (S);
  2006. length (<< rand(2,2); rand(10,10); 3.14 >>);
  2007.  
  2008. printf ("\tlength");
  2009.  
  2010. log (A);
  2011. log ([A]);
  2012. log (3.14);
  2013.  
  2014. log (Z);
  2015. log ([Z]);
  2016. log (3.14j);
  2017.  
  2018. printf ("\tlog");
  2019.  
  2020. log10 (A);
  2021. log10 ([A]);
  2022. log10 (3.14);
  2023.  
  2024. log10 (Z);
  2025. log10 ([Z]);
  2026. log10 (3.14j);
  2027.  
  2028. printf ("\tlog10");
  2029.  
  2030. max (A);
  2031. max ([A]);
  2032. max (3.14);
  2033.  
  2034. max (Z);
  2035. max ([Z]);
  2036. max (3.14j);
  2037.  
  2038. printf ("\tmax");
  2039.  
  2040. maxi (A);
  2041. maxi ([A]);
  2042. maxi (3.14);
  2043.  
  2044. maxi (Z);
  2045. maxi ([Z]);
  2046. maxi (3.14j);
  2047.  
  2048. printf ("\tmaxi");
  2049. printf ("\n");
  2050.  
  2051. members (<< rand(2,2); rand(); rand(3,3) >>);
  2052.  
  2053. printf ("\tmembers");
  2054.  
  2055. min (A);
  2056. min ([A]);
  2057. min (3.14);
  2058.  
  2059. min (Z);
  2060. min ([Z]);
  2061. min (3.14j);
  2062.  
  2063. printf ("\tmin");
  2064.  
  2065. mini (A);
  2066. mini ([A]);
  2067. mini (3.14);
  2068.  
  2069. mini (Z);
  2070. mini ([Z]);
  2071. mini (3.14j);
  2072.  
  2073. printf ("\tmini");
  2074.  
  2075. mod (A,A);
  2076. mod ([A],A);
  2077. mod (3.14,2);
  2078.  
  2079. mod (Z,Z);
  2080. mod ([Z],Z);
  2081. mod (3.14j,2);
  2082.  
  2083. printf ("\tmod");
  2084.  
  2085. nan ();
  2086.  
  2087. ones(3,3);
  2088. ones([4,4]);
  2089.  
  2090. printf ("\tones");
  2091.  
  2092. prod (A);
  2093. prod ([A]);
  2094. prod (3.14);
  2095.  
  2096. prod (Z);
  2097. prod ([Z]);
  2098. prod (3.14j);
  2099.  
  2100. printf ("\tprod");
  2101.  
  2102. real (A);
  2103. real ([A]);
  2104. real (3.14);
  2105.  
  2106. real (Z);
  2107. real ([Z]);
  2108. real (3.14j);
  2109.  
  2110. printf ("\treal");
  2111. printf ("\n");
  2112.  
  2113. reshape (A, A.nr*A.nc, 1);
  2114. reshape (A, [A.nr*A.nc], [1]);
  2115. reshape (Z, Z.nr*Z.nc, 1);
  2116. reshape (Z, [Z.nr*Z.nc], [1]);
  2117.  
  2118. printf ("\treshape");
  2119.  
  2120. round (A);
  2121. round ([A]);
  2122. round (3.14);
  2123.  
  2124. round (Z);
  2125. round ([Z]);
  2126. round (3.14j);
  2127.  
  2128. printf ("\tround");
  2129.  
  2130. show(A);
  2131. show(Z);
  2132. show([A]);
  2133. show([Z]);
  2134.  
  2135. printf ("\tshow");
  2136.  
  2137. sign (A);
  2138. sign ([A]);
  2139. sign (3.14);
  2140.  
  2141. sign (Z);
  2142. sign ([Z]);
  2143. sign (3.14j);
  2144.  
  2145. printf ("\tsign");
  2146.  
  2147. sin (A);
  2148. sin ([A]);
  2149. sin (3.14);
  2150.  
  2151. sin (Z);
  2152. sin ([Z]);
  2153. sin (3.14j);
  2154.  
  2155. printf ("\tsin");
  2156.  
  2157. size (A);
  2158. size ([A]);
  2159. size (3.14);
  2160.  
  2161. size (Z);
  2162. size ([Z]);
  2163. size (3.14j);
  2164.  
  2165. printf ("\tsize");
  2166.  
  2167. sizeof (A);
  2168. sizeof ([A]);
  2169. sizeof (3.14);
  2170.  
  2171. sizeof (Z);
  2172. sizeof ([Z]);
  2173. sizeof (3.14j);
  2174.  
  2175. printf ("\tsizeof");
  2176. printf ("\n");
  2177.  
  2178. sort (A);
  2179. sort ([A]);
  2180. sort (3.14);
  2181.  
  2182. sort (Z);
  2183. sort ([Z]);
  2184. sort (3.14j);
  2185.  
  2186. printf ("\tsort");
  2187.  
  2188. sqrt (A);
  2189. sqrt ([A]);
  2190. sqrt (3.14);
  2191.  
  2192. sqrt (Z);
  2193. sqrt ([Z]);
  2194. sqrt (3.14j);
  2195.  
  2196. printf ("\tsqrt");
  2197.  
  2198. srand ();
  2199. srand ("clock");
  2200. srand (SEED);
  2201.  
  2202. printf ("\tsrand");
  2203.  
  2204. strsplt (S);
  2205. printf ("\tstrsplt");
  2206.  
  2207. strtod ("string");
  2208. strtod (S);
  2209. strtod ("1.2");
  2210. strtod (["1.2", "1e3"]);
  2211.  
  2212. printf ("\tstrtod");
  2213.  
  2214. sum (A);
  2215. sum ([A]);
  2216. sum (3.14);
  2217.  
  2218. sum (Z);
  2219. sum ([Z]);
  2220. sum (3.14j);
  2221.  
  2222. printf ("\tsum");
  2223.  
  2224. tan (A);
  2225. tan ([A]);
  2226. tan (3.14);
  2227.  
  2228. tan (Z);
  2229. tan ([Z]);
  2230. tan (3.14j);
  2231.  
  2232. printf ("\ttan");
  2233. printf ("\n");
  2234.  
  2235. tmpnam ();
  2236. printf ("\ttmpnam");
  2237.  
  2238. type (A);
  2239. type ([A]);
  2240. type (Z);
  2241. type ([Z]);
  2242. type (3.14);
  2243. type (S);
  2244. type ("string");
  2245.  
  2246. printf ("\ttype");
  2247.  
  2248. zeros (3,3);
  2249. zeros ([3], [3]);
  2250. zeros ([3,3]);
  2251.  
  2252. printf ("\tzeros");
  2253.  
  2254. printf ("\n");
  2255.  
  2256. srand(SEED);
  2257. rand("default");
  2258.  
  2259. //
  2260. //-------------------------- Test rand () -----------------------------
  2261. //
  2262. rand ("normal", 5, 1);
  2263. xrand = rand(4000,1);
  2264.  
  2265. mean = function(x)
  2266. {
  2267.   local(m);
  2268.  
  2269.   m = size (x)[1];
  2270.   if( m == 1 ) 
  2271.   { 
  2272.     m = size (x)[2];
  2273.   }
  2274.  
  2275.   return sum( x ) / m;
  2276. };
  2277.  
  2278. std = function(x)
  2279. {
  2280.   local(i, m, s);
  2281.  
  2282.   if(class(x) != "num") { error("std() requires NUMERICAL input"); }
  2283.  
  2284.   m = x.nr;
  2285.   if( m == 1 ) 
  2286.   { 
  2287.     return sqrt( sum( (x - mean(x)) .^ 2 ) / (x.nc - 1) );
  2288.   else
  2289.     for( i in 1:x.nc) {
  2290.       s[i] = sqrt( sum( (x[;i] - mean(x[;i])) .^ 2 ) / (x.nr - 1) );
  2291.     }
  2292.     return s;
  2293.   }
  2294. };
  2295.  
  2296. if (!(mean (xrand) > 4.9 && mean (xrand) < 5.1)) 
  2297.   { error ("error in random"); }
  2298. if (!(std (xrand) > 0.9 && std (xrand) < 1.1))
  2299.   { error ("error in random"); }
  2300. printf("\tpassed rand test...\n");
  2301.  
  2302. rand("default");
  2303.  
  2304. //
  2305. //-------------------------- Test norm () -----------------------------
  2306. //
  2307.  
  2308. tn = [1,2,3,4;2,1,2,3;3,2,1,2;4,3,2,1  ];
  2309. if (norm(tn,"m") != 4 ) { error ("incorrect norm computation"); }
  2310. if (norm(tn,"1") != 10) { error ("incorrect norm computation"); }
  2311. if (norm(tn,"i") != 10) { error ("incorrect norm computation"); }
  2312. printf("\tpassed norm test...\n");
  2313.  
  2314. //
  2315. //-------------------------- Test qr () ------------------------------
  2316. //
  2317.  
  2318. a = ohess(4);
  2319. qa = qr (a);
  2320. if (max (max (abs (qa.q*qa.r - a)))/(X*norm (a)*a.nr) > eps)
  2321.   { error ("possible qr() problems"); }
  2322.  
  2323. z = ohess (4) + ohess(4)*1i;
  2324. qz = qr (z);
  2325. if (max (max (abs (qz.q*qz.r -  z)))/(X*norm (z)*z.nr) > eps)
  2326.   { error ("possible qr() problems"); }
  2327.  
  2328. printf("\tpassed qr test...\n");
  2329.  
  2330. //
  2331. //-------------------------- Test schur () ----------------------------
  2332. //
  2333.  
  2334. a = randsvd (10, 10);
  2335. sa = schur (a);
  2336. if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
  2337.   { error ("possible schur() problems"); }
  2338.  
  2339. z = rand (4,4) + rand(4,4)*1i;
  2340. sz = schur (z);
  2341. if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
  2342.   { error ("possible schur() problems"); }
  2343.  
  2344. a = randsvd (10, -10);
  2345. sa = schur (a);
  2346. if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
  2347.   { error ("possible schur() problems"); }
  2348.  
  2349. z = rand (4,4) + rand(4,4)*1i;
  2350. sz = schur (z);
  2351. if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
  2352.   { error ("possible schur() problems"); }
  2353.  
  2354. printf("\tpassed schur test...\n");
  2355.  
  2356. //
  2357. //-------------------------- Test schord () ----------------------------
  2358. //
  2359.  
  2360. s2a = schord (sa, 2, 4);
  2361. if (max (max (abs (s2a.z*s2a.t*s2a.z' - a)))/(X*norm (a)*a.nr) > eps)
  2362. { 
  2363.   error ("possible schord() problems"); 
  2364. }
  2365.  
  2366. s2z = schord (sz, 3, 1);
  2367. if (max (max (abs (s2z.z*s2z.t*s2z.z' - z)))/(X*norm (z)*z.nr) > eps)
  2368. { 
  2369.   error ("possible schord() problems"); 
  2370. }
  2371.  
  2372. printf("\tpassed schord test...\n");
  2373.  
  2374. //
  2375. //-------------------------- Test chol () -----------------------------
  2376. //
  2377.  
  2378. c = lehmer(10);
  2379. u = chol (c);
  2380. if (max (max (abs (u'*u - c)))/(X*norm (c)*c.nr) > eps)
  2381. { 
  2382.   error ("possible chol() problems"); 
  2383. }
  2384.  
  2385. cz = lehmer(10) + lehmer(10)*1j;
  2386. cz = symm(cz);
  2387. uz = chol (cz);
  2388. if (max (max (abs (uz'*uz - cz)))/(X*norm (cz)*cz.nr) > eps)
  2389. { 
  2390.   error ("possible chol() problems"); 
  2391. }
  2392.  
  2393. printf("\tpassed chol test...\n");
  2394.  
  2395. //
  2396. //--------------------------- Test inv () --------------------------------
  2397. //
  2398.  
  2399. a = randsvd(10,10);
  2400. b = ones(10,1);
  2401. x = inv(a) * b;
  2402. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  2403. { 
  2404.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2405.   printf ("\tA*X - B:\n");
  2406.   abs (a*x - b)
  2407.   error ("possible inv() problems\n");
  2408. }
  2409.  
  2410. az = randsvd(10,10) + randsvd(10,10)*1j;
  2411. bz = rand(10,1) + rand(10,1)*1j;
  2412. xz = inv (az)*bz;
  2413. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2414. { 
  2415.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2416.   printf ("\tA*X - B:\n");
  2417.   abs (az*xz - bz)
  2418.   error ("possible inv() problems\n");
  2419. }
  2420.  
  2421. printf("\tpassed inv test...\n");
  2422.  
  2423. //
  2424. //-------------------------- Test solve () -----------------------------
  2425. //
  2426.  
  2427. //
  2428. // Real - General case
  2429. //
  2430.  
  2431. a = randsvd(10,10);
  2432. b = ones(10,1);
  2433. x = solve (a,b);
  2434. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  2435. { 
  2436.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2437.   printf ("\tA*X - B:\n");
  2438.   abs (a*x - b)
  2439.   error ("possible solve() problems\n");
  2440. }
  2441.  
  2442. //
  2443. // Real - Symmetric case
  2444. //
  2445.  
  2446. s = symm (randsvd(10,10));
  2447. b = ones(10,1);
  2448. x = solve (s,b);
  2449. if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
  2450. { 
  2451.   printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  2452.   printf ("\tA*X - B:\n");
  2453.   abs (s*x - b)
  2454.   error ("possible solve() problems\n");
  2455. }
  2456.  
  2457. //
  2458. // Complex - General  case
  2459. //
  2460.  
  2461. az = randsvd(10,10) + randsvd(10,10)*1j;
  2462. bz = rand(10,1) + rand(10,1)*1j;
  2463. xz = solve (az,bz);
  2464. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2465. { 
  2466.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2467.   printf ("\tA*X - B:\n");
  2468.   abs (az*xz - bz)
  2469.   error ("possible solve() problems\n");
  2470. }
  2471.  
  2472. //
  2473. // Complex - Symmetric case
  2474. //
  2475.  
  2476. sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
  2477. bz = rand(10,1) + rand(10,1)*1j;
  2478. xz = solve (sz,bz);
  2479. if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
  2480. { 
  2481.   printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  2482.   printf ("\tA*X - B:\n");
  2483.   abs (sz*xz - bz)
  2484.   error ("possible solve() problems\n");
  2485. }
  2486.  
  2487. printf("\tpassed solve test...\n");
  2488.  
  2489. //
  2490. //-------------------------- Test factor() / backsub() -----------------------------
  2491. //
  2492.  
  2493. //
  2494. // Real - General case
  2495. //
  2496.  
  2497. a = randsvd(10,10);
  2498. b = ones(10,1);
  2499. f = factor (a);
  2500. x = backsub (f,b);
  2501. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  2502. { 
  2503.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2504.   printf ("\tA*X - B:\n");
  2505.   abs (a*x - b)
  2506.   error ("possible factor/backsub problems\n");
  2507. }
  2508.  
  2509. //
  2510. // Real - Symmetric case
  2511. //
  2512.  
  2513. s = symm (randsvd(10,10));
  2514. f = factor (s);
  2515. x = backsub (f,b);
  2516. if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
  2517. { 
  2518.   printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  2519.   printf ("\tA*X - B:\n");
  2520.   abs (s*x - b)
  2521.   error ("possible factor/backsub problems\n");
  2522. }
  2523.  
  2524. s = symm (randsvd(10,10));
  2525. f = factor (s, "s");
  2526. x = backsub (f,b);
  2527. if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
  2528. { 
  2529.   printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  2530.   printf ("\tA*X - B:\n");
  2531.   abs (s*x - b)
  2532.   error ("possible factor/backsub problems\n");
  2533. }
  2534.  
  2535. //
  2536. // Complex - General  case
  2537. //
  2538.  
  2539. az = randsvd(10,10) + randsvd(10,10)*1j;
  2540. bz = rand(10,1) + rand(10,1)*1j;
  2541. fz = factor (az);
  2542. xz = backsub (fz,bz);
  2543. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2544. { 
  2545.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2546.   printf ("\tA*X - B:\n");
  2547.   abs (az*xz - bz)
  2548.   error ("possible factor/backsub problems\n");
  2549. }
  2550.  
  2551. az = randsvd(10,10) + randsvd(10,10)*1j;
  2552. bz = rand(10,1) + rand(10,1)*1j;
  2553. fz = factor (az, "g");
  2554. xz = backsub (fz,bz);
  2555. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2556. { 
  2557.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2558.   printf ("\tA*X - B:\n");
  2559.   abs (az*xz - bz)
  2560.   error ("possible factor/backsub problems\n");
  2561. }
  2562.  
  2563. //
  2564. // Complex - Symmetric case
  2565. //
  2566.  
  2567. sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
  2568. fz = factor(sz);
  2569. xz = backsub (fz,bz);
  2570. if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
  2571. { 
  2572.   printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  2573.   printf ("\tA*X - B:\n");
  2574.   abs (sz*xz - bz)
  2575.   error ("possible factor/backsub problems\n");
  2576. }
  2577.  
  2578. sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
  2579. fz = factor(sz, "s");
  2580. xz = backsub (fz,bz);
  2581. if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
  2582. { 
  2583.   printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  2584.   printf ("\tA*X - B:\n");
  2585.   abs (sz*xz - bz)
  2586.   error ("possible factor/backsub problems\n");
  2587. }
  2588.  
  2589. printf("\tpassed factor/backsub test...\n");
  2590.  
  2591. //
  2592. //------------------------------ Test lu() ---------------------------------
  2593. //
  2594.  
  2595. static (swap);
  2596.  
  2597. lu = function ( A )
  2598. {
  2599.   local (i, l, u, pvt, x)
  2600.  
  2601.   if (A.nr != A.nc) { error ("lu() requires square A"); }
  2602.  
  2603.   x = factor (A, "g");    // Do the factorization
  2604.  
  2605.   //
  2606.   // Now create l, u, and pvt from lu and pvt.
  2607.   //
  2608.  
  2609.   l = tril (x.lu, -1) + eye (size (x.lu));
  2610.   u = triu (x.lu);
  2611.   pvt = eye (size (x.lu));
  2612.  
  2613.   //
  2614.   // Now re-arange the columns of pvt
  2615.   //
  2616.  
  2617.   for (i in 1:max (size (x.lu)))
  2618.   {
  2619.     pvt = pvt[ ; swap (1:pvt.nc, i, x.pvt[i]) ];
  2620.   }
  2621.   return << l = l; u = u; pvt = pvt >>;
  2622. };
  2623.  
  2624. //
  2625. //  In vector V, swap elements I, J
  2626. //
  2627.  
  2628. swap = function ( V, I, J )
  2629. {
  2630.   local (v, tmp);
  2631.   v = V;
  2632.   tmp = v[I];
  2633.   v[I] = v[J];
  2634.   v[J] = tmp;
  2635.   return v;
  2636. };
  2637.  
  2638. a = randsvd(10,10);
  2639. lua = lu (a);
  2640. if (max (max (abs(a - lua.pvt*lua.l*lua.u)))/(X*norm (a)*a.nr) > eps)
  2641. { 
  2642.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2643.   printf ("\tA - p*l*u:\n");
  2644.   abs (a - lua.pvt*lua.l*lua.u)
  2645.   error ("possible lu()/factor() problems\n");
  2646. }
  2647.  
  2648. //
  2649. // Real
  2650. az = randsvd(10,10) + randsvd(10,10)*1j;
  2651. luz = lu (az);
  2652. if (max (max (abs (az - luz.pvt*luz.l*luz.u)))/(X*norm (az)*az.nr) > eps)
  2653. { 
  2654.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2655.   printf ("\tA - p*l*u:\n");
  2656.   abs (az - luz.ovt*luz.l*luz.u)
  2657.   error ("possible lu()/factor()() problems\n");
  2658. }
  2659.  
  2660. printf("\tpassed lu/factor test...\n");
  2661.  
  2662. //
  2663. //-------------------------- Test svd ()   -----------------------------
  2664. //
  2665.  
  2666. a = randsvd(10,10);
  2667. s = svd (a);
  2668. if (max (max (abs (s.u*diag(s.sigma)*s.vt - a)))/(X*norm (a)*a.nr) > eps)
  2669. {
  2670.   error ("possible svd() problems"); 
  2671. }
  2672.  
  2673. z = randsvd(10,10) + rand(10,10)*1j;
  2674. sz = svd (z);
  2675. if (max (max (abs (sz.u*diag(sz.sigma)*sz.vt - z)))/(X*norm (z)*z.nr) > eps)
  2676. { 
  2677.   error ("possible svd() problems"); 
  2678. }
  2679.  
  2680. printf("\tpassed svd test...\n");
  2681.  
  2682. //
  2683. //-------------------------- Test hess ()  -----------------------------
  2684. //
  2685.  
  2686. a = randsvd(10,10);
  2687. h = hess (a);
  2688. if (max (max (abs (h.p*h.h*h.p' - a)))/(X*norm (a)*a.nr) > eps)
  2689. { 
  2690.   error ("possible hess() problems");
  2691. }
  2692.  
  2693. z = randsvd(10,10) + randsvd(10,10)*1j;
  2694. hz = hess (z);
  2695. if (max (max (abs (hz.p*hz.h*hz.p' - z)))/(X*norm (z)*z.nr) > eps)
  2696. { 
  2697.   error ("possible hess() problems"); 
  2698. }
  2699.  
  2700. printf("\tpassed hess test...\n");
  2701.  
  2702. //
  2703. //-------------------------- Test lyap () ------------------------------
  2704. //
  2705.  
  2706. lyap = function ( A, B, C )
  2707. {
  2708.   local (A, B, C)
  2709.  
  2710.   if (!exist (B)) 
  2711.   { 
  2712.     B = A';    // Solve the special form: A*X + X*A' = -C
  2713.   }
  2714.  
  2715.   if ((A.nr != A.nc) || (B.nr != B.nc) || (C.nr != A.nr) || (C.nc != B.nr)) {
  2716.     error ("Dimensions do not agree.");
  2717.   }
  2718.  
  2719.   //
  2720.   // Schur decomposition on A and B
  2721.   //
  2722.  
  2723.   sa = schur (A);
  2724.   sb = schur (B);
  2725.  
  2726.   //
  2727.   // transform C
  2728.   //
  2729.  
  2730.   tc = sa.z' * C * sb.z;
  2731.  
  2732.   X = sylv (sa.t, sb.t, tc);
  2733.  
  2734.   //
  2735.   // Undo the transformation
  2736.   //
  2737.  
  2738.   X = sa.z * X * sb.z';
  2739.  
  2740.   return X;
  2741. };
  2742.  
  2743. a = randsvd (10,10);
  2744. b = rand (10,10);
  2745. c = rand (10,10);
  2746.  
  2747. x = lyap (a, b, c);
  2748. if (max (max (abs (a*x + x*b + c)))/(X*norm(c)*norm(a)*norm(b)) > eps)
  2749. { 
  2750.   error ("possible problems with lyap() or sylv()"); 
  2751. }
  2752.  
  2753. printf("\tpassed lyap test...\n");
  2754.  
  2755. //
  2756. //-------------------------- Test eig () ------------------------------
  2757. //
  2758.  
  2759. trace = function(m) 
  2760. {
  2761.   local(i, tr);
  2762.  
  2763.   if(m.class != "num") { 
  2764.     error("must provide NUMERICAL input to trace()");
  2765.   }
  2766.  
  2767.   tr = 0;
  2768.   for(i in 1:min( [m.nr, m.nc] )) {
  2769.     tr = tr + m[i;i];
  2770.   }
  2771.  
  2772.   return tr;
  2773. };
  2774.  
  2775. eye = function( m , n ) 
  2776. {
  2777.   local(i, N, new);
  2778.  
  2779.   if (!exist (n))
  2780.   {
  2781.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  2782.     new = zeros (m[1], m[2]);
  2783.     N = min ([m[1], m[2]]);
  2784.   else
  2785.     if (class (m) == "string" || class (n) == "string") {
  2786.       error ("eye(), string arguments not allowed");
  2787.     }
  2788.     if (max (size (m)) == 1 && max (size (n)) == 1)
  2789.     {
  2790.       new = zeros (m[1], n[1]);
  2791.       N = min ([m[1], n[1]]);
  2792.     else
  2793.       error ("matrix arguments to eye() must be 1x1");
  2794.     }
  2795.   }
  2796.   for(i in 1:N)
  2797.   {
  2798.     new[i;i] = 1.0;
  2799.   }
  2800.   return new;
  2801. };
  2802.  
  2803. //
  2804. // Standard eigenvalue problem
  2805. //
  2806.  
  2807. a = randsvd(10,10);
  2808. ta = trace (a);
  2809. sa = symm (a);
  2810. tsa = trace (sa);
  2811.  
  2812. z = randsvd(10,10) + randsvd(10,10)*1i;
  2813. tz = trace (z);
  2814. sz = symm (z);
  2815. tsz = trace (sz);
  2816.  
  2817. tol = 1.e-6;
  2818.  
  2819. if (!(ta < sum(eig(a).val) + tol && ta > sum(eig(a).val) - tol))
  2820. {
  2821.   error ("error in eig");
  2822. }
  2823. if (!(tsa < sum(eig(sa).val) + tol && tsa > sum(eig(sa).val) - tol))
  2824. {
  2825.   error ("error in eig");
  2826. }
  2827. if (abs(tz)+tol < abs(sum(eig(z).val)) && abs(tz)+tol > abs(sum(eig(z).val)))
  2828. {
  2829.   error ("error in eig");
  2830. }
  2831. if (abs(tsz)+tol < abs(sum(eig(sz).val)) && abs(tsz) > abs(sum(eig(sz).val)))
  2832. {
  2833.   error ("error in eig");
  2834. }
  2835.  
  2836. //
  2837. // Generalized eigenvalue problem
  2838. //
  2839.  
  2840. b = randsvd(10,10);
  2841. sb = symm (b) + eye(size(b))*3;
  2842. tb = trace (b);
  2843. tsb = trace (sb);
  2844.  
  2845. zb = randsvd(10,10) + randsvd(10,10)*1i;
  2846. szb = symm (zb) + eye(size(zb))*3;
  2847. tzb = trace (zb);
  2848. tszb = trace (szb);
  2849.  
  2850. eig(a,b);    // not sure of a good way to check these yet
  2851. eigs(sa,sb);
  2852. eigs(sa,sb);
  2853.  
  2854. eig(z, zb);
  2855. eigs(sz, szb);
  2856. eigs(sz, szb);
  2857.  
  2858. printf("\tpassed eig test...\n");
  2859.  
  2860. //
  2861. //-------------------------- Test fft () -----------------------------
  2862. //
  2863.  
  2864. if (100 != fft(ones(100,1))[1]) { error ("error in fft()"); }
  2865. printf("\tpassed fft test...\n");
  2866.  
  2867. //
  2868. //------------------------- Fibonacci Test -------------------------
  2869. //
  2870. //  Calculate Fibonacci numbers
  2871. //
  2872.  
  2873. i=1; 
  2874. while ( i < 2 ) { 
  2875.   i=i+1;
  2876.   a=0; b=1;
  2877.   while ( b < 10000 ) {
  2878.     c = b;
  2879.     b = a+b;
  2880.     a = c;
  2881.   }
  2882. }
  2883. if ( b != 10946 ) {
  2884.   error("failed fibonacci test");
  2885. else
  2886.   printf("\tpassed fibonacci test...\n");
  2887. }
  2888.  
  2889. //
  2890. //------------------------- Factorial Test -------------------------
  2891. //
  2892.  
  2893. fac = function(a) 
  2894. {
  2895.   if(a <= 1) 
  2896.   {
  2897.     return 1
  2898.   else
  2899.     return a*$self(a-1)
  2900.   }
  2901. };
  2902.  
  2903. if(fac(10) != 3628800)
  2904. { 
  2905.   error(); else printf("\tpassed factorial test...\n"); 
  2906. }
  2907.  
  2908. //
  2909. //--------------------------- ACK Test ----------------------------
  2910. //
  2911.  
  2912. ack = function(a, b) 
  2913. {
  2914.   if(a == 0) { return b + 1; }
  2915.   if(b == 0) { return $self(a-1, 1); }
  2916.   return $self(a-1, $self(a, b-1));
  2917. };
  2918.  
  2919. if(ack(2,2) != 7) 
  2920. {
  2921.   error("error in ack() test");
  2922.   else
  2923.   printf("\tpassed ACK test...\n");
  2924. }
  2925.  
  2926. //
  2927. //------------------------- Prime Test -----------------------------
  2928. //
  2929. // An example that finds all primes less than limit
  2930. //
  2931.  
  2932. primes = function (limit) 
  2933. {
  2934.   local(prime, cnt, i, j, k);
  2935.   
  2936.   i = 1; cnt = 0;
  2937.   for(k in 2:limit) 
  2938.   {
  2939.     j = 2;
  2940.     while(mod(k,j) != 0) 
  2941.     {
  2942.       j++;
  2943.     }
  2944.     if(j == k)             // Found prime
  2945.     {
  2946.       cnt++;
  2947.       prime[i;1] = k;
  2948.       i++;
  2949.     }
  2950.   }
  2951.   return prime;
  2952. };
  2953.  
  2954. if(max(size(primes(100))) != 25) 
  2955. {  
  2956.   error("error in prime test");
  2957.   else
  2958.   printf("\tpassed prime test...\n");
  2959. }
  2960.  
  2961. //
  2962. //--------------------------- Fib Min Test -----------------------------
  2963. //
  2964. //    fibmin() will minimize an arbitrary function 
  2965. //    in 1D using Fibonacci search
  2966.  
  2967. f065 = function(x)
  2968. {
  2969.   return (x - 0.65) * (x - 0.65);
  2970. };
  2971.  
  2972. fib = function(x)
  2973. {
  2974.   local (n, a, b);
  2975.   
  2976.   a = 1;
  2977.   b = 1;
  2978.   if (x >= 2)
  2979.   {
  2980.     n = x - 1;
  2981.     for (n in n:1:-1)
  2982.     {
  2983.       c = b;
  2984.       b = a + b;
  2985.       a = c;
  2986.       n = n - 1;
  2987.     }
  2988.   }
  2989.   return b;
  2990. };
  2991.  
  2992. //  Minimize a 1D function using Fibonacci search
  2993. //  f = function to minimize
  2994. //  xlo = lower bound
  2995. //  xhi = upper bound
  2996. //  n = number of iterations (the bigger the more accurate)
  2997.  
  2998. fibmin = function(f, xlo, xhi, n) 
  2999. {
  3000.   local(a, b, x, y, ex, ey, k, lo, hi);
  3001.   
  3002.   lo = xlo;
  3003.   hi = xhi;
  3004.   k = n;
  3005.   for (k in k:2:-1)
  3006.   {
  3007.     a = fib(k - 2) / fib(k);
  3008.     b = fib(k - 1) / fib(k);
  3009.     x = lo + (hi - lo) * a;
  3010.     y = lo + (hi - lo) * b;
  3011.     ex = f(x);
  3012.     ey = f(y);
  3013.     if (ex >= ey)
  3014.     {
  3015.       lo = x;
  3016.       else
  3017.       hi = y;
  3018.     }
  3019.     //  printf("%d: (%g %g) %g %g %g %g\n",  k, a, b, lo, hi, ex, ey);
  3020.   }
  3021.   return (lo + hi) / 2;
  3022. };
  3023.  
  3024. //
  3025. // Simple example using above function to mimize f065. Answer is 0.65
  3026. //
  3027.  
  3028. x = fibmin(f065, 0, 1, 30); // printf("f(%g)=%g\n", x, f065(x));
  3029. if (abs(x - 0.65) > 1e-6)
  3030. {
  3031.   printf("x = %f\n", x);
  3032.   error("failed fibmin test");
  3033. }
  3034.  
  3035. printf("\tpassed fibmin test...\n");
  3036.  
  3037. //
  3038. //--------------------- Nasty Function Test ------------------------
  3039. //
  3040.  
  3041. printf("\tStarting Nasty Function Test...");
  3042. printf("\tthis will take awhile\n");
  3043. check = function( a, b, c, d, e, f, g, h ) 
  3044. {
  3045.   if ( a+b+c+d == e+f+g+h && ...
  3046.        a^2+b^2+c^2+d^2 == e^2+f^2+g^2+h^2 && ...
  3047.        a^3+b^3+c^3+d^3 == e^3+f^3+g^3+h^3 ) 
  3048.   {
  3049.     return 1;
  3050.   else
  3051.     return 0;
  3052.   }
  3053. };
  3054.  
  3055. for(a in 8:10) {
  3056.   for(b in 7:(a-1)) {
  3057.     for(c in 6:(b-1)) {
  3058.       for(d in 5:(c-1)) {
  3059.         for(e in 4:(d-1)) {
  3060.           for(f in 3:(e-1)) {
  3061.             for(g in 2:(f-1)) {
  3062.               for(h in 1:(g-1)) {                  
  3063.               if(check( a, b, c, d,  e, f, g, h ) || ...
  3064.                      check( a, e, c, d,  b, f, g, h ) || ...
  3065.                      check( a, f, c, d,  e, b, g, h ) || ...
  3066.                      check( a, g, c, d,  e, f, b, h ) || ...
  3067.                      check( a, h, c, d,  e, f, g, b ) || ...
  3068.                      check( a, b, e, d,  c, f, g, h ) || ...
  3069.                      check( a, b, f, d,  e, c, g, h ) || ...
  3070.                      check( a, b, g, d,  e, f, c, h ) || ...
  3071.                      check( a, b, h, d,  e, f, g, c ) || ...
  3072.                      check( a, b, c, e,  d, f, g, h ) || ...
  3073.                      check( a, b, c, f,  e, d, g, h ) || ...
  3074.                      check( a, b, c, g,  e, f, d, h ) || ...
  3075.                      check( a, b, c, h,  e, f, g, d ) || ...
  3076.                      check( a, e, f, d,  b, c, g, h ) || ...
  3077.                      check( a, e, g, d,  b, f, c, h ) || ...
  3078.                      check( a, e, h, d,  b, f, g, c ) || ...
  3079.                      check( a, f, g, d,  e, b, c, h ) || ...
  3080.                      check( a, f, h, d,  e, b, g, c ) || ...
  3081.                      check( a, g, h, d,  e, f, b, c ) || ...
  3082.                      check( a, b, e, f,  c, d, g, h ) || ...
  3083.                      check( a, b, e, g,  c, f, d, h ) || ...
  3084.                      check( a, b, e, h,  c, f, g, d ) || ...
  3085.                      check( a, b, f, g,  e, c, d, h ) || ...
  3086.                      check( a, b, f, h,  e, c, g, d ) || ...
  3087.                      check( a, b, g, h,  e, f, c, d ) || ...
  3088.                      check( a, e, f, g,  e, f, g, h ) || ...
  3089.                      check( a, e, f, h,  e, f, g, h ) || ...
  3090.                      check( a, e, g, h,  e, f, g, h ) || ...
  3091.                      check( a, f, g, h,  e, f, g, h ) ) { cnt++; }
  3092.               }
  3093.             }
  3094.           }
  3095.         }
  3096.       }
  3097.     }
  3098.   }
  3099. }
  3100.  
  3101. if(1) {  // figure out the value of cnt, and check!
  3102.   printf("\tpassed nasty function test...\n");
  3103. else
  3104.   error();
  3105. }
  3106.  
  3107. //
  3108. //------------------ Test More Advanced Functions --------------------
  3109. //
  3110.  
  3111. printf( "\tStarting the lqr/ode test..." );
  3112. printf( "\tthis will take awhile\n" );
  3113.  
  3114. lqr = function( a, b, q, r ) 
  3115. {
  3116.   local( k, s,...
  3117.          m, n, mb, nb, mq, nq,...
  3118.          e, v, d );
  3119.  
  3120.   m = size(a)[1]; n = size(a)[2];
  3121.   mb = size(b)[1]; nb = size(b)[2];
  3122.   mq = size(q)[1]; nq = size(q)[2];
  3123.   
  3124.   if ( m != mq || n != nq ) 
  3125.   {
  3126.     fprintf( "stderr", "A and Q must be the same size.\n" );
  3127.     quit
  3128.   }
  3129.     
  3130.   mr = size(r)[1]; nr = size(r)[2];
  3131.   if ( mr != nr || nb != mr ) 
  3132.   {
  3133.     fprintf( "stderr", "B and R must be consistent.\n" );
  3134.     quit
  3135.   }
  3136.     
  3137.   nn = zeros( m, nb );
  3138.     
  3139.   // Start eigenvector decomposition by finding eigenvectors of Hamiltonian:
  3140.     
  3141.   e = eig( [ a, solve(r',b')'*b'; q, -a' ] );
  3142.   v = e.vec; d = e.val;
  3143.     
  3144.   index = sort( real( d ) ).ind;
  3145.   d = real( d[ index ] );
  3146.   
  3147.   if ( !( d[n] < 0 && d[n+1] > 0 ) ) 
  3148.   {
  3149.     fprintf( "stderr", "Can't order eigenvalues.\n" );
  3150.     quit
  3151.   }
  3152.   
  3153.   chi = v[ 1:n; index[1:n] ];
  3154.   lambda = v[ (n+1):(2*n); index[1:n] ];
  3155.   s = -real(solve(chi',lambda')');
  3156.   k = solve( r, nn'+b'*s );
  3157.   
  3158.   return << k=k; s=s >>;
  3159.   
  3160. };
  3161.  
  3162. // Now run a little test problem.
  3163.  
  3164. k = 1; m = 1; c = .1;
  3165. a = [0     ,1    ,0    , 0;
  3166.     -k/m, -c/m,  k/m,  c/m;
  3167.      0,     0,    0,    1;
  3168.      k/m,  c/m, -k/m, -c/m ];
  3169. b = [ 0; 1/m; 0; 0 ];
  3170. qxx = diag( [0, 0, 100, 0] );
  3171. ruu = [1];
  3172. K = lqr( a, b, qxx, ruu ).k;
  3173.  
  3174. dot = function( t, x ) 
  3175. {
  3176.   global (a, b, K)
  3177.   return (a-b*K)*x + b*K*([1,0,1,0]');
  3178. };
  3179.  
  3180. x = ode ( dot, 0, 15, [0,0,0,0], .02, 1e-5, 1e-5 );
  3181.  
  3182. m = maxi( x[;2] );
  3183.  
  3184. if ( (abs( x[m;2] - 1.195 ) > 0.001)  || ...
  3185.      any (abs( x[x.nr;2,4] - 1 ) > 0.001) ) 
  3186. {
  3187.   printf( "\tfailed***\n" );
  3188.   else
  3189.   printf( "\tpassed the lqr/ode test...\n" );
  3190. }
  3191.  
  3192. printf("Elapsed time = %10.3f seconds\n", toc() );
  3193. "FINISHED TESTS"
  3194.